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ABSTRACT 


l  nmanned  Air  \ehicles  have  become  increasingly  important  on  the  modern 
battlefield.  The  restrictive  requirement  for  runways  and  special  equipment  to  take  off 
and  land  was  partially  solved  by  the  vertical  take  off  and  landing  Airborne  Remotely 
Operated  Device,  AROD.  Work  done  at  the  Naval  Postgraduate  School  has  modified 
the  AROD  to  not  only  land  and  launch  vertically,  but  to  fly  horizontally  for  the 
majority  of  the  mission.  To  realize  these  capabilities,  as  well  as  that  of  autonomous 
flight,  an  accurate  computer  model  was  required  of  both  the  AROD  and  the  avionics 
test  bed  aircraft.  Bluebird,  in  order  to  design  the  control  and  navigation  systems.  High 
fidelity,  non-linear  equations  of  motion  were  derived  in  matrix  form  that  represented 
any  six  degree  of  freedom  aircraft  model,  and  were  then  tailored  for  use  on  specific 
aircraft.  Computer  modeling  of  the  resulting  equations  of  motion,  as  well  as  the 
sensors  used  on  the  aircraft,  was  done  using  SlMt  LINK  and  M.-VTLAB  software.  The 
resulting  computer  model  provides  a  non  -linear  system  of  equations,  which  are  easily 
linearized  at  any  desired  flight  condition,  as  required  by  the  proposed  control  and 
navigation  system  design. 
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I.  INTRODUCTION 


A.  IMPORTANCE  OF  UNMANNED  AIR  VEHICLES 

Unmanned  Aerial  Vehicles  have  become  increasingly  more  important,  both  on 
the  battlefield  and  in  civilian  service,  since  the  Ryan  Q-2C  Firebee  target  drone 
introduced  the  “modern  age”  of  UAYs  in  1960  [Ref.  Siu  91].  From  that  time  on, 
military  planners  have  assured  that  UAVs  have  the  capability  to  collect  intelligence, 
target  er  my  positions,  gather  bomb  damage  assessment,  as  well  as  perform  many 
other  tasks.  The  real  benefit  in  using  unmanned  aircraft  lies  in  the  fact  that  m  ny 
missions  can  be  performed  deep  in  enemy  territory,  all  without  endangering  the  lives 
of  pilots,  or  risking  the  loss  of  a  much  more  expensive  aircraft.  With  the  recent  use  of 
UAVs  in  Operation  Desert  Storm,  improvements  in  the  current  technology  are  both 
indicated  and  desirable. 

The  most  capable  UAV  in  service  of  the  United  States  Navy  and  Marine  Corps 
today  is  the  i  ioneer  Short  Range  UAV.  The  system  is  hampered,  though,  by  the 
large  amount  of  runway  and  special  equipment  needed  to  launch  and  land  the  aircraft. 
These  requirements  limit  the  usefulness  of  the  Pioneer  by  keeping  the  aircraft  take¬ 
off  and  landing  area  well  away  from  the  areas  where  the  ground  forces  are  operating. 
This  distance  then  leads  to  longer  transit  times  to  and  from  the  assigned  operating 
area,  and  thus  a  shorter  time  on  station.  However,  what  is  really  required  in  many 
instances  by  the  g.ound  forces  is  an  aircraft  that  can  respond  quickly  to  a  changing 
tactical  situation.  The  Airborne  Remotely  Operated  Device,  AROD,  was  an  attempt 
by  Sandia  National  Laboratories  [Ref.  Wh  87],  in  response  to  a  requirement  by  ihe 
Naval  Ocean  Systems  Center.  NOSC,  to  respond  to  these  needs. 


The  United  States  Marine  Corps  had  set  a  requirement  for  a  short  range,  direct, 
support  UAV  as  described  in  [Ref.  MCG  87]: 

“  ...to  allow  the  front  line  commander  to  see  “over  the  next  hill",  to  a 

distance  of  two  kilometers  ...” 

The  AROD  was  designed  to  be  a  ducted  fan,  hovering  device  carrying  a  fiber  optic 
data  link  and  on  board  cameras.  AROD  testing  was  canceled  [Ref.  Sa  89]  as  the 
Department  of  Defense  requirements  grew,  requiring  a  minimum  range  of  30  km  for 
all  Short  Range  UAVs.  The  AROD  was  incapable  of  this  kind  of  range,  however:  the 
design  is  still  potentially  useful. 

The  Unmanned  Air  Vehicle  Flight  Research  Lab,  UAV  FRL,  at  the  Naval  Post¬ 
graduate  School  has  proposed  a  solution  using  the  AROD  that  would  satisfy  the 
DoD  short  range  UAV  requirements,  while  maintaining  the  important  capability  for 
vertical  take-offs  and  landings. 

B.  UAV  RESEARCH  USING  THE  AROD  VEHICLE 

Unmanned  Aerial  Vehicle  research  underway  at  NPS  has  taken  the  AROD  air¬ 
frame  and  fitted  it  with  wings  from  the  Aquila  UAV  [Ref.  Kre  92.  Sto  93]  in  order 
to  give  the  AROD  forward  flight  capability.  The  proposed  configuration  will  give  the 
Archytas  (an  AROD  with  wings)  the  ability  to  take-off  and  land  vertically  and  then 
transition  to  horizontal  flight  for  the  mission.  This  design  will  explore  new  technol¬ 
ogy,  driven  by  the  goals  established  by  the  UAV  Joint  Project  Office  [Ref.  DOD  92] 
of: 

•  Take  off  weight  under  200  lb 

•  Carry  a  50  lb  payload 


•  Fly  at  a  maximum  speed  of  150  kts 


•  Take-off  and  land  in  an  area  30m  by  60m 

The  vertical  flight  is  accomplished  with  a  powerful  ducted  fan. which  causes  a  great 
deal  of  gyroscopic  coupling  and  torque  when  producing  enough  thrust  to  lift  the 
aircraft.  Therefore,  in  order  to  achieve  stable  take-offs  and  landings,  a  three-axis 
autopilot  is  a  necessary  feature  of  the  aircraft.  Additional  capabilities  desired  in 
the  final  version  of  the  Archvtas  are  guidance  and  navigation  systems  which  will 
allow  autonomous  operation,  as  well  as  a  Global  Positioning  System  aided  autoland 
capability. 

C.  REQUIREMENT  FOR  MODELING 

Simulation  and  modeling  of  the  aircraft  are  essential  to  the  successful  design 
of  a  control  system  capable  of  autonomous  flight.  The  model  must  be  a  very  high 
fidelity,  non-linear  model,  that  can  be  easily  linearized  at  any  given  flight  condition. 
The  model  should  be  able  to  interpolate  between  data  points  resulting  from  wind 
tunnel  testing  in  order  to  simulate  the  highly  non-linear  transition  from  vertical  to 
horizontal  flight.  Moreover,  the  model  must  also  be  capable  of  including  the  outputs 
of  the  sensors  as  inputs  to  the  control  and  navigation  system  for  sensors  located  at 
any  arbitrary  location  on  the  aircraft. 

This  thesis  develops  a  six  degree  of  freedom  model  for  the  AROD  in  the  vertical 
flight  regime,  as  well  as  for  an  aircraft  in  a  fixed  wing  configuration.  This  test  aircraft, 
named  Bluebird,  is  used  to  test  the  guidance,  navigation,  and  control,  GNC,  systems 
in  horizontal  flight,  since  there  currently  is  no  aerodynamic  data  available  for  the 
Archytas  configuration.  Use  of  the  Bluebird  will  provide  the  capability  to  design  and 
test  a  GNC  system  on  a  stable  aircraft  before  the  first  flight  on  the  Archytas. 
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II.  BACKGROUND 


A.  DESCRIPTION  OF  AROD 

The  AROD  was  designed  by  the  Sandia  Research  Laboratory  in  Albuquerque, 
New  Mexico  in  a  project  managed  by  NOSC.  The  vehicle  possessed  no  flying  surfaces 
and  relied  solely  on  powered  lift  for  flight.  Control  of  the  aircraft  was  obtained  through 
the  use  of  four  fixed  anti-torque  vanes  and  four  moveable  control  vanes  positioned  in 
the  propeller  wash  of  the  duct  [Ref.  We  88].  The  main  features  of  the  AROD  were 
vertical  take-off  and  landing.  VTOL.  flight,  lightweight  construction,  compact  size, 
and  minimal  support  equipment  required.  However,  the  AROD  required  most  of  the 
engine  output  to  maintain  the  powered  lift,  so  very  little  excess  thrust  was  left  for 
translational  flight. 

An  important  aspect  of  the  AROD  design  was  the  improvement  in  static  perfor¬ 
mance  provided  by  the  efficiency  of  the  ducted  fan  design.  The  addition  of  the  shroud 
around  the  three-bladed  propeller  resulted  i.i  increased  mass  flow  through  the  fan, 
and  thus  more  static  thrust  when  compared  to  a  conventional  propeller  configuration 
[Ref.  Kre  92].  The  AROD  is  shown  in  Figure  2.1  and  characteristics  of  the  AROD 
are  tabulated  in  Table  2.1.  The  moveable  control  vanes  are  all  used  in  combination 
to  exert  the  desired  control  forces  on  AROD.  Roll  control  is  achieved  by  deflecting  the 
four  vanes  in  the  same  direction,  while  pitch  and  yaw  control  is  obtained  by  deflecting 
a  pair  of  vanes  in  the  required  direction.  The  numbering  of  the  vanes  is  shown  in 
Figure  2.2  and  the  combinations  of  vane  deflections  required  for  positive  roll,  pitch, 
and  yaw  motion  are  given  in  Table  2.2. 
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Figure  2.1:  Airborne  Remotely  Operated  Device,  [Ref.  Siu  91] 

B.  AROD  MODELING 

The  AROD  vehicle  has  been  the  subject  of  several  theses  at  NPS.  The  designs 
presented  in  the  theses  rely  on  the  AROD  model  given  by  Sandia  Labs  in  the  original 
design  of  their  controller  [Ref.  Wh  87,  Wh  91].  This  model  was  based  on  the  more 
classical  technique  of  linearizing  an  aircraft  model,  based  on  dimensional  derivatives 
in  a  state  space  form.  The  resulting  model  was  acceptable  for  the  AROD  in  a  hovering 
and  near  vertical  translational  flight  mode,  but  was  not  easily  adaptable  to  anything 
other  than  the  narrow  range  of  conditions  planned  for  AROD.  The  Sandia  Labs 
papers  also  pointed  out  several  types  of  coupling  in  the  AROD.  The  most  prominent 
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TABLE  2.1:  PHYSICAL  CHARACTERISTICS  OF  AROD 


Inlet  Diameter,  A 

29.25  in 

Propeller  Radius,  R 

12  in 

Exit  Radius 

23.375  in 

Inlet  Area  Ratio 

1.219 

Exit  Area  Ratio 

1.115 

Exterior  Contour 

Tapered  Rear 

Propeller  Location,  %  chord 

25  % 

Number  of  Blades 

3 

Engine  Speed,  Max. 

8000  rpm 

Engine  Speed,  Nom. 

6500  rpm 

Tip  Speed,  Max. 

838  fpm 

Tip  Speed,  Nom. 

680  fpm 

7.25  HP/ f2 

Mass  Moment  of  Inertia.  IT 

Mass  Moment  of  Inertia.  /„ 

mumi 

Mass  Moment  of  Inertia.  L 

Prop  Mass  Moment  of  Inertia.  lTjr 

li 

Prop  Mass  Moment  of  Inertia,  Iry 

Prop  Mass  Moment  of  Inertia.  ITZ 

WWWMtnEHMl 

TABLE  2.2:  VANE  DEFLECTION  COMBINATIONS  FOR  POSITIVE 
ANGLES 


Vane  Combination 

Roll,  $ 

V,  +  V2  +  V3  4  V4 

Pitch.  0 

\  2  -  v4 

Yaw,  $ 

Vi  -  v  3 

of  the  coupling  effects  is  the  gyroscopic  coupling  between  the  pitch  and  yaw  axes 
resulting  from  the  large  amount  of  angular  momentum  contributed  to  the  aircraft 
by  the  propeller.  Another  dynamic  coupling  exists  between  the  altitude-rate  and  the 
vehicle  attitude,  since  a  loss  of  lift  due  to  thrust  will  occur  when  the  vehicle  is  tilted  to 
generate  horizontal  motion.  Yet  a  third  dynamic  coupling  exists  between  the  altitude 
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Figure  2.2:  AROD  Direction  of  Positive  Vane  Deflections 


and  roll  control  loops,  since  the  reactive  torques  applied  to  the  roll  axis  vary  as  the 
engine  speed  is  varied.  Sandia  Labs  also  provided  data  for  modeling  both  the  engine 
and  the  servos  as  second  order  transfer  functions  which  were  used  in  this  thesis. 

Additional  information  was  obtained  by  Weir  [Ref.  We  88]  in  wind  tunnel  test¬ 
ing.  This  information  included  non-dimensional  derivatives  for  vane  effectiveness  and 
non-dimensional  stability  derivatives.  The  report  also  stated  that  the  control-vane 
effectiveness  is  constant  out  to  at  least  25  deg  of  deflection.  Wind  tunnel  data  were 
also  presented  to  show  that  control  vane  effectiveness  in  approximately  the  same  for 
translational  flight  as  for  hovering  flight.  This  equivalence  is  due  to  the  fact  that  the 
vanes  are  located  in  the  high  speed  flow  aft  of  the  propeller  and  are  not  significantly 
affected  by  the  freestream. 


TABLE  2.3:  PHYSICAL  CHARACTERISTICS  OF  BLUEBIRD 


Weight 

55  lbs 

Average  Wing  Chord,  c 

1.802  / 

WTing  Span,  b 

12.42  / 

Planform  Surface  Area,  S 

Engine  Power 

4.0  HP 

Mas:  Moment  of  Inertia,  IT 

Mass  Moment  of  Inertia.  Iy 

KSBSI 

Mass  Moment  of  Inertia,  lz 

wmmam 

C.  DESCRIPTION  OF  BLUEBIRD 

The  Bluebird  aircraft  was  acquired  as  a  test  bed  for  guidance  and  navigation 
systems.  Ultimately,  these  systems  will  be  installed  on  the  Archvtas  aircraft.  The 
Bluebird  is  a  conventional  aircraft  that  will  be  used  to  test  systems  for  a  similar 
configuration  to  the  Archytas  in  forward  flight.  The  aircraft  model  was  developed 
in  the  same  manner  as  for  AROD,  as  will  be  described  in  Chapter  III.  Physical 
characteristics  for  the  Bluebird  are  given  in  Table  2.3. 
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III.  AIRCRAFT  EQUATIONS  OF  MOTION 


A.  NOTATION 

In  this  section  the  notation  used  in  modeling  the  equations  of  motion  is  intro¬ 
duced.  This  notation  is  common  in  the  field  of  robotics  (see  [Ref.  Sil  91]  and  [Ref. 
Cra  86]).  and  is  shown  below  in  Figure  3.1.  The  following  conventions  are  used: 


Figure  3.1:  Relative  Position  of  Coordinate  Systems 


•  {A}  represents  the  coordinate  system  with  basis  vectors.  xa-Vai  and  z^. 

•  aPq  represents  the  position  of  point  Q.  expressed  in  {A}. 

•  aVq  represents  the  velocity  of  point  Q,  measured  in  {.4}  and  expressed  in  {.4}. 


B[AVq)  represents  the  velocity  of  point  Q.  measured  in  {.4},  and  expressed  ii 

im. 


•  gR  is  a  rotation  matrix,  also  called  a  direction  cosine  matrix.  A  free  vector 
in  one  coordinate  system,  that  is  a  vector  that  can  be  mov(d  parallel  to  itself 
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unthout  change  e.g.,  bVq  can  be  expressed  in  another  coordinate  system  by 
using  the  rotation  matrix: 

A(BvQ)  =  imBvQ) 

•  aSIb  is  the  angular  velocity  of  the  { B }  coordinate  system  with  respect  to  {-4}, 
and  expressed  in  {A}. 

■  •  b(aQb)  is  the  angular  velocity  of  { B },  with  respect  to  {A}.  and  expressed  in 

{B}- 

•  Additional  information  can  be  added  to  the  subscripts  i.e.,  A  Pgo  is  the  position 
of  the  origin  of  {/?},  expressed  in  {A}. 

B.  COORDINATE  SYSTEMS 

In  order  to  derive  equations  of  motion  for  a  rigid  airplane,  the  following  coor¬ 
dinate  systems  and  assumptions  will  be  used: 

•  {U}  represents  the  inertial  tangent  plane  coordinate  system  attached  to  Earth. 

•  {£?}  represents  the  body  fixed  coordinate  system. 

•  All  sensors  are  located  at  the  c.g.  (This  assumption  will  be  lifted 
in  a  later  section) 

•  The  mass  of  the  aircraft  remains  constant. 

•  Given  a  vector  v,  its  derivative  with  respect  to  {B}  is  denoted  as  jt(-) 
and 

its  derivative  with  respect  to  {U}  is  denoted  as  (') 
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The  {B}  coordinate  system  is  defined  with  A g  as  the  thrust  axis.  A  positive  roll 
rate,  p.  is  clockwise  when  looking  in  the  positive  X  direction.  The  positive  direction 
for  Eg,  the  pitch  axis,  is  out  the  right  wing  .  A  positive  pitch  rate,  q,  is  defined  as 
clockwise  looking  in  the  positive  V  direction.  The  Zg  axis  is  the  yaw  axis,  and  a 
positive  yaw  rate,  r,  is  defined  as  clockwise,  looking  in  the  positive  Z  direction. 

To  simplify  the  notation  in  places  where  it  becomes  cumbersome.  The  following 
definitions  are  introduced: 

•  vq  represents  the  velocity  of  an  arbitrary  point.  Q,  measured  and  expressed  in 

{O- 

•  *'sc  represents  the  velocity  of  the  origin  of  {/?}.  measured  and  expressed  in 
{l  },  i.e.  1  VBO  =  t’fio- 

•  bvq  represents  the  velocity  of  point  Q.  measured  in  {['}  and  expressed  in  {B}. 
i.e.  V’lb)  =  %. 

•  *t,'g  represents  the  angular  velocity  of  {£?}.  measured  and  expressed  in  {l  ’}.  i.e. 
1  ftg  =  u :g. 

•  Bujg  represents  the  angular  velocity  of  {B}  measured  in  {l  },  and  expressed  in 
{B}.  i.e.  e(f  nB)  = 

C.  SPATIAL  ORIENTATION 
1.  Euler  Angles 

The  spatial  orientation  of  a  rigid  body  [Ref.  Ju  92]  can  be  defined  by  the 
three  Euler  angles,  $.0,  and  called  roll,  pitch  and  yaw  and  defined  in  Figure  3.2. 
The  Euler  angles, in  turn,  can  be  used  to  define  a  rotation  between  two  coordinate 
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Figure  3.2:  Z-Y-X  Euler  Angle  Rotation  Sequence 

systems.  This  rotation  is  obtained  using  Euler's  theorem: 

Any  number  of  rotations  about  different  axes  through  a  point  must,  in 
the  end,  remain  equivalent  to  a  single  rotation. 

For  the  case  of  conventional  aircraft,  a  3-2-1  rotation  sequence  is  used  [Ref.  Sch  92], 
where  the  aircraft  is  yawed,  pitched,  then  rolled.  In  the  cases  investigated  here,  0  is 
small,  and  in  steady  state  flight  is  equal  to  the  angle  of  attack,  a.  $  can  be  expected 
to  be  anywhere  from  ±  60deg  in  normal  flight  and  can  be  anywhere  from  ±  180deg  in 
acrobatic  flight.  represents  the  heading  of  the  aircraft  and  of  course  can  range  from 
0  to  360  deg.  The  transformation  from  inertial  coordinates{f to  body  coordinates 
{5},  is  carried  out  as  follows,  and  is  shown  in  Figure  3.2. 

1.  The  inertial  coordinate  system  is  represented  by  the  vector  L  V',  with  the  com¬ 
ponents  x.  y,  and  z.  The  first  rotation  is  made  about  the  z  axis  through  an  angle 
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ty.  Now  the  vector  is  expressed  as  2V  with  the  components  x2.y2.  and  z2-  Since 
the  rotation  was  about  the  z  axis,  the  z2  component  remains  unchanged.  The 
resulting  elemental  matrix  is: 

cos  'I*  sin  0 

AZ('P)  =  —  sin  $  cos 'I'  0  (3.1) 

0  0  1 

2.  Now  the  rotation  is  made  about  the  new  y  axis,  y2,  through  an  angle  0.  This 
results  in  a  third  coordinate  system  with  the  vector  expressed  as  3V  ,  and  having 
components  x3,y3,  and  z3.  This  rotation  does  not  change  the  y3  component. 

The  resulting  elemental  matrix  is: 

cos  0  0  —  sin  0 

A/(0)  =010.  (3.2) 

sin0  0  cos  0 

3.  Lastly,  the  rotation  about  the  x3  axis  through  an  angle  4>  is  made  to  give  the 

vector  expressed  in  body  coordinates.  B\'.  Now  the  resulting  matrix  is 

'1  0  0 

M ($)  =  0  cos$  sin  $  (3.3) 

0  —  sin  4>  cos<l> 


When  the  matrices  are  multiplied  together  in  the  correct  sequence.  A/(<J>)A/(0)A/(  ^ ). 
the  result  is  the  BH  direction  cosine  matrix,  expressed  in  terms  of  Euler  angles  as 
shown 

cos'I,cos0  sin^cosO  —  sin0 

cos'L  sin0  sin3>  —  si  n't  cos<l>  sin0  sin$  sin'L  +  costy  cos^>  cosOsinO  .  (3.4) 

cos^  sin0  cos4>  +  sin'I' sin<I>  sin0  cos$  sin^  —  costy  sin<J>  cos0cos<l> 

The  next  step  is  to  develop  the  kinematic  differential  equations  that  de¬ 
scribe  the  change  in  Euler  angles  with  time.  Following  the  method  used  in  [Ref. 
Sch  92],  the  matrix  of  differential  equations,  ft,  can  be  written  as  a  sum  of  individual 


Euler  angle  rates: 


ft  =  A/(4>) 


+  A/(4>)A/(0) 


+  A/  ( <I> )  A/  ( 0 )  A/  ( ) 


(3.5) 


\ 

When  the  matrix  algebra  in  Equation  3.5  is  done,  the  resulting  kinematic  differential 
equations  for  p.q. and  r  are  given  as: 

p  1  0  —  sin  0  4> 

q  =  0  cos  4>  cos©  sin  4>  0  (3.6) 

r  0  —  sin  cos  0  cos 

The  matrix  on  the  right  hand  side  of  Equation  3.6  is  invertible  for  all  0  ^  t/2,  and 

can  be  used  to  solve  for  the  Euler  angle  rates.  0  and 

4>  1  sin  $  tan©  cos  $  tan  0  p 

0  =  0  cos  $  —  sin  4>  q  t3.7) 

0  sin  $  sec  0  cos4>  sec  0  r 

Bv  integrating  Equation  3.7,  the  time  history  of  the  Euler  angles  can  be  obtained. 

The  Euler  angle  method  has  one  drawback.  In  the  kinematic  differential 

equations  derived,  a  singularity  occurs  for  some  particular  value  of  either  <J>.  0.  or  ty. 

In  Equation  3.".  this  singularity  occurs  at  0  =  ±tt/2.  which  means  ihe  that  the 

coordinate  transformation  is  useful  for  an  aircraft  in  horizontal  flight,  but  useless 

for  an  aircraft  which  requires  a  vertical  take-off  or  extended  periods  of  vertical  flight. 

Thus,  another  type  of  transformation  is  necessary  for  aircraft  that  spend  considerable 

time  flying  at  conditions  near  0  =  tt/2. 

The  first  alternative  is  simply  to  use  another  one  of  the  12  possible  Euler 

angle  transformations.  It  has  been  shown  that  the  rotation  matrix.  R.  is  made  up 

of  sequential  rotations  and  can  be  characterized  as  the  product  o'-  three  individual 

matrices  where 

R  =  Mj(02)  A/o(0i).  (3.8) 

where  the  rotation  sequence  (a,  0,  q)  repiesents  one  of  the  combinations  cf  integers 

(1.2.3) ,  (1.3.2),  (1,2,1).  (1.3.1) 

(2. 1.3) .  (2, 3.1).  (2. 1.2),  (2.  3,  2) 

(3.2. 1).  (3. 1.2).  (3. 2. 3).  (3, 1.3) 
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and  the  M,{Oj)  are  the  rotation  matrices 


r  i 

0 

0 

AM*)  = 

0  cos  9  sin 

9 

.0  - 

sin 

9  cos 

9  _ 

cos  9 

0 

—  sin 

9  ' 

AM*)  = 

0 

1 

0 

sin  9 

0 

cos  9 

cos  9 

:in  9 

0  ' 

AM*)  = 

—  sin 

9 

cos  9 

0 

0 

0 

1 

One  can  see  that  by  substituting  in  (3,2,1),  and  tf.O,  and  $  for  9\,92,  and  03  re¬ 
spectively,  the  matrix  shown  in  Equation  3.4  is  obtained.  Euler  angle  transformation 
matrices  for  each  of  these  combinations  have  been  calculated  and  are  tabulated  in 


[Ref.  Ka  83].  The  best  Euler  angle  rotation  sequence  for  an  aircraft  with  flight  at 

0  =  r/2  was  selected  as  a  2-3-1,  or  Y-Z-X,  sequence.  This  sequence  will  allow  flight 

at  0  =  7r/2  with  no  corresponding  singularity  in  the  kinematic  differential  equations. 

Note  that  this  matrix  is  invertible  for  ^  7r/2.  The  2-3-1  rotation  matrix.  R,  is 

described  in  terms  of  Euler  angles  as 

cos  0  cos 'I'  sin 'I'  —  sin  0  cos 

—  cos  0  sin  cos  *I>  4-  sin  4>  sin  0  cos  cos  $  sin  0  sin  T  cos  $  +  sin  0  cos  0 
cos  0  sin  'P  sin  $  +  cos  sin  9  —  cos  'P  sin  0  —  sin  0  sin  4n  <!>  -f  cos  0  cos  $ 

(3.10) 


The  kinematic  differential  equations  can  be  found  in  [Ref.  Ka  83]  as 


'  <D  ' 

0 

1 

cos  $ 

.  V  . 

1 

0 

0 


—  cos  4>  sin  ^ 
cos  $ 

sin  $  cos  $ 


sin  4>  sin  'P 
—  sin 
cos  $  cos 


P 

Q 

r 


(3.11) 


These  equations  can  now  be  integrated  to  find  the  time  history  of  the  Euler  angles. 


2.  Quaternions 


Another  choice  for  the  expression  of  a  body's  spatial  orientation  is  the 
use  of  quaternions.  Quaternions  eliminate  the  disadvantage  of  the  singularity  in  the 
second  rotation  that  is  associated  with  the  Euler  angles.  Quaternions  have  been  in 


15 


use  for  quite  some  time,  having  been  discovered  by  Euler  in  a  search  for  complex 
numbers  [Ref.  Mo  84).  A  quaternion  is  defined  as  [Ref.  Mo  84]: 


9  —  To  +  i  Ti  +  j $2  +  k  83 , 


(3.12) 


where  the  parameters  are  repr  sented  by  various  authors  as  S,  a,  b,  c  by  [Ref.  Ro  58], 

\,  £.  t).  and  C  by  [Ref.  Whi  59],  and  q4,  <?,.  q2,  and  q3  by  [Ref.  Sil  91]. 

The  components  To,  /?i ,  /?2-  and  T3  are  real  numbers  and  the  terms  i,j,  and  k 

are  defined  in  the  typical  manner  for  complex  numbers,  where 

i2  =  —  1  ij  =  -ji  =  k 
J2  =  “I  jk  =  ~kj  =  ? 
k2  =  —  1  ki  =  —ik  =  j 

The  norm  of  a  quaternion,  q’q ,  is  required  to  be  1: 

9  9*  =  q'  q  =  8l  +  8\  +  =  1 , 

since  9*  =  T0  -  ?T1  -  j0 2  -  A-T3. 

It  can  be  shown  that  R  can  be  represented  as  follows: 


(3.13) 


where 


= 

’  rll  r12  ^13 

r21  r22  r23 

r  31  r32  r33 

Hi  = 

■ 

-  T2  - 
>?2 

r12  = 

2(  Tl  02  + 

^0  T3  ) 

r13  = 

2(T,T3- 

8082) 

r21  = 

2(T./?2  - 

8083) 

r22  = 

81-01  +  0]  - 

r22  =  -  T,  +  T|  -  T32  .  (3.14) 

r23  =  2(T2T3  +  To$l) 

r3i  =  2(Ti  T3  4-  T0T2) 

r32  =  ‘2(02  03  ~  8o8\) 

r33  =  /?o  -  0?  ~  T|  +  8% 

All  that  is  required  now  is  to  determine  the  kinematic  differential  equations  using  the 
quaternions. 
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By  expressing  subsequent  rotations  from  one  coordinate  system  to  another, 
where  3'  orients  Fi  to  F,  and  (3"  orients  F2  to  Fj,  an  algebraic  approach  [Ref.  Mo  8-1] 
can  be  used  to  relate  F2  to  F. 


R(3)  =  R(0")R(0') 

RtAS)  =  EL,  R,k(3")Rkj(3')  ' 


(3.15) 

Now  the  0's  can  be  expressed  in  terms  of  /3'  s  and  3"' s  with  the  following  result 

0=R(0")0\  (3.16) 

where 


R(3)  = 


(3.17) 


$0  —ft\  —32 

3\  3o  3$  —  32 

32  —3$  3o  3\ 

L  &  32  -3x  00  J 

By  regarding  the  second  rotation  in  Equation  3.16as  infinitesimal,  the  following  result 


is  obtained 


where 


3  =  -R(u:*)3, 


(3.18) 


3  = 


3o 

01 

02 

k 


and 


R{  w*)  = 

Equation  3.18  can  be  rewritten  as 


3o 

3  = 

3\ 

02 

and 

LJ*  = 

.  03  . 

'  0 

—  u.’l 

— u-’2 

—^3 

U>1 

0 

u^3 

—  a,’2 

— u>3 

0 

uq 

. 

U>2 

-U.’, 

0  . 

0 

P 

<7 

r 


(3.19) 


(3.20) 


0=-R(0)<*. 


(3.21) 


and  in  this  form  can  be  integrated  to  give  the  time  history  of  the  orientation  of  the 
body. 
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Using  quaternions  has  the  following  advantages  over  Euler  angles  in  repre¬ 
senting  spatial  orientation  of  a  rigid  body: 

•  Four  states  required  to  express  the  spatial  orientation. 

•  Requires  almost  30  9c  fewer  calculations  [Ref.  Ro  58].  mainly  because  no  non¬ 
linear,  trigonometric  equations  need  to  be  calculated. 

•  No  singularities  in  Equation  3.21  at  any  body  attitude. 

D.  DERIVATION  OF  EQUATIONS  OF  MOTION 

The  derivation  of  equations  of  motion  fora  general  six  degree  of  freedom  airplane 
model  can  be  divided  into  two  parts.  The  first  part  is  simply  the  determination 
of  the  equations  of  motion  for  any  rigid  body  in  space.  It  is  dependent  only  on 
the  linear  and  angular  momenta  of  the  body.  The  second  part  is  the  calculation 
of  aerodynamic,  gravitational,  and  thrust  forces  on  the  airplane.  These  forces  are 
particular  to  a  certain  aircraft  and  in  general  can  be  represented  by  the  stability  and 
control  derivatives  described  later  in  the  thesis. 
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1.  Linear  Equations 

Equations  for  linear  motion  can  be  calculated  using  Newton's  law.  F  —  m  a. 
Since  most  of  the  sensor  information  available  for  feedback  to  the  control  and  nav¬ 
igation  systems  is  available  in  the  {B}  coordinate  system,  the  terms  for  linear  ac¬ 
celerations,  as  well  as  forces  and  moments,  will  be  expressed  in  the  body  coordinate 
system.  First  the  position  of  the  aircraft  c.g.  is  determined  as  1  Pbo-  Then  Coriolis’ 
theorem  is  applied  to  obtain  linear  velocities  for  the  aircraft.  Coriolis'  theorem  is 
then  reapplied  to  derive  the  expression  for  linear  accelerations.  Then 

1  y bo  =  1  Pbo  (3.22) 

Both  sides  of  Equation  3.22  are  prcmultiplied  by  B  R  to  get: 

BRL  Vbo  —  BRl  Pbo 

or 

Bvbo  =  BPbo  (3.23) 

Now  consider  Coriolis'  theorem 

d 

A  =  —A+*'xA,  (3.24) 

dt 

where  A  and  jtA  use  the  notation  for  derivatives  previously  defined  in  Section  A.  The 
term  u:  x  ,4  represents  the  difference  between  the  relative  velocity  as  measured  from 
the  rotating  and  non-rotating  axes  [Ref.  (Jre  88]. 

Equation  3.24  is  applied  to  bvbo  in  Equation  3.23  to  get: 

Bi'B°  =  -jjB<-'BO  +  B^'b  x  Bvbo ■  (3.25) 

Newton’s  law  can  now  be  written  as 

1  F  =  m  1  a 

=  rn  vbo •  (3.26) 
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s 

where  1  F  is  the  total  external  force  applied  to  the  aircraft  c.g.  Equation  3.26  is 
premuhiplied  by  BR  to  obtain  the  result: 

t  —  m  vn  i'bo 

=  rn  Bi'BO •  (3.27) 

when  Bi’Bo  is  substituted  into  Equation  3.27,  the  final  result  for  B F  is 

B  r  m  ^  i  B  ,  v-  B  \ 

t  =  m  ( —  I'BO  +  x  vBO ) 
at 

=  m  —bi'bo  +  m  Bu.’b  x  Bi'so-  (3.28) 

at 

2.  Angular  Equations 

The  equations  for  angular  accelerations  are  derived  using  Euler's  law  for 
preservation  of  angular  momentum.  These  equations  are  also  derived  for  the  aircraft 
c.g.  by  applying  Coriolis'  theorem  to  the  equation  for  Euler's  law: 

VLbo  =  vSBO •  (3.29) 

where  1  Lbo  is  the  angular  momentum  of  the  aircraft  and  L  Sbo  is  the  total  ex¬ 
ternal  moment  applied  to  the  aircraft  c.g.  Euler's  law  can  be  rewritten  in  { B }  by 
premultiplying  Equation  3.29  by  BR  to  get 

B LBo  -  f-R1  Xbo-  (3.30) 

Using  Coriolis'  theorem  in  Equation  3.24.  B Lbo  can  be  rewritten  as 

B Lbo  =  ~rB Lbo  +  b^b  *  B Lbo-  (3.31) 

at 

The  angular  momentum,  B Lbo,  is  defined  as  the  product  of  the  inertia  tensor,  1b, 
and  the  body's  angular  velocity,  b^b,  or 

S/— 7  B  ,  if  B  , 

L  —  Ib  w’b  +  IB  u,r. 


20 


where  Ir  and  Bu)r  are  the  moment  of  inertia  and  the  angular  velocity  of  the  propeller, 
respectively.  When  this  term  is  substituted  into  Equation  3.31,  the  result  is 

B Lbo  =  -t.{1bB^b  +  1rB*>r)  +  Bu*  x  [1  bB^'b  +  IrB^'r)-  (3.32) 

at 

where  7g  is  the  inertia  tensor  for  the  aircraft  and  Jr  is  the  inertia  tensor  for  any 
significant  spinning  object  on  the  aircraft,  such  as  a  propeller,  turbine  ,  or  other 
rotating  disk.  The  term,  Bu;g,  is  the  angular  velocity  of  the  rotor,  expressed  in  {5}. 
We  can  carry  out  the  differentiation  in  Equation  3.32  to  get 

B Lbo  =  1b-tB^b  +  Irj-b^'r  +  B^b  x  {IbB^b  +  Irb^>r)-  (3.33) 
at  at 

However,  since  d/dt(Bu.’g)  =  bJu'b  and  d/dt(Bu,'R)  is  assumed  to  be  very  small,  Equa¬ 
tion  3.33  results  in 

B Lbo  =  IbB(L>b  +  Bu -b  x  (7sBu.’g  +  IrB^'r)  (3.34) 

Now  the  result  in  Equation  3.34  can  be  substituted  into  Equation  3.29: 

B Nbo  =  IbB<L>b  +  x  (IbBub  +  1rB^'r)-  (3.35) 

The  term  IrBur  can  be  disregarded  if  it  is  insignificant  compared  to  7g  and  Bu>B 
[Ref.  Ros  79]. 

3.  State  Equations 

In  the  preceeding  sections,  kinematic  equations  for  the  motion  of  a  rigid 
body  were  derived  in  matrix  form.  These  equations  can  be  assembled  into  a  state 
space  representation  of  the  kinematic  equations  of  motion.  First,  Equations  3.28  and 
3.35  can  be  written  as 


Equation  3.36  can  be  rearranged  to  yield 


d 

rn  B  i'bo 

m  {—b^'b  x  Bi'bo ) 

+bF)  ' 

It 

Ibb^'b 

—  Bu.'b  *  {Ibb^’b  +  1rb*-'r) 

+B.Y)  _ 

The  terms  on  the  left  hand  side  of  Equation  3.37  can  be  normalized  by  multiplying 
by  1/m  and  /g1.  with  the  final  result  of 


d 

Bi'bo 

-b^b  x  bvbo 

+  ¥  ’ 

dt 

— /g1  Bu.'b  x  (Ibb^b  +  1rb^r) 

+  ip  ®-vJ 

(3.38) 


E.  EXTERNAL  FORCES  AND  MOMENTS 


Section  D.  gives  the  equations  for  kinematic  motion,  as  shown  in  Equation  3.38, 
for  any  rigid  body.  Now  it  is  necessary  to  distinguish  between  the  different  platforms 
to  be  modeled  in  order  to  give  an  accurate  representation  of  the  aircraft.  This  is 
achieved  by  computing  the  actual  BF  and  B X  acting  on  the  aircraft.  These  forces 
and  moments  are  those  due  to  gravitational,  propulsive,  and  aerodynamic  effects, 
written  as 


‘  BF  ‘ 

B.V 

BFgRAV  +  B  ^  PROP  +  Bt.AERO 
B ^PROP  +  B  A  aero 


(3.39) 


1.  Aerodynamic  Forces  and  Moments 

The  aerodynamic  force  and  moment  terms  are  determined  by  using  a  first 
order  Taylor  series  expansion  around  a  given  nominal  operating  point.  This  operating 
point  is  the  state  chosen  to  represent  the  aircraft's  flight  condition.  Each  term  in  the 
series  is  a  partial  derivative  of  B F  or  B X  with  respect  to  the  aerodynamic  variables. 
u/l\  q.  3 ,  p.  q.  r  [Ref.  Sch  92.  Th  89]: 


Faero  =  6FT‘jc'  +  SFyx'  +  +  F0.  (3.40) 


Similarly,  moment  terms  can  be  written  as 


A aero  —  F\z‘.r'  +  +  FN^A  +  Aq. 


(3.41) 


where  x'  is  given  by 


and  x'  is  given  as 


T-  =  r”  3  a  3L  _±i 

[V'I  '2U'2V'2IJI 


x'  =  [0,a]. 

Control  inputs  are  represented  by  the  vector  A: 

A  =  [6',6r,6a] 

where  6e,  6r ,  and  6a  are  the  elevator,  rudder,  and  aileron  inputs,  respectively. 
Equations  3.40-3.44  can  now  be  combined  as  follows: 


_e(dC  ,80.,  80  „  . 

qS^dx'T  +  8i'X  +5AA  +  Cfo}’ 


(3.42) 


(3.43) 


(3.44) 


(3.45) 


where  q  =  l/2pl/2,  5  =  diag {5, 5, 5, 56, 5c,  56},  and  C  is  the  matrix  of  non- 
dimensional  stability  derivatives  differentiated  with  respect  to  the  terms  defined  in 
Equation  3.42,  3.43,  or  3.44.  §p  is  defined  as: 


cLv 

Cl, 

cLa 

Clp 

Cl , 

CL, 

CYu 

Cyb 

cYa 

Cyp 

C-Yn 

cYr 

Cdv 

Cd , 

Cdo 

CoP 

Cd, , 

Cdt 

c,v 

Ci, 

C,a 

CiP 

Ci, 

Clr 

c 

^  my 

Cm , 

Cma 

Cmp 

r 

^  mq 

Cm, 

r 

nil 

r 

'-rip 

Cna 

CnP 

Cnq 

Cm, 

§p  is  very  similar  to  |p,  except  that  only  the  q  and  [3  terms  are  normally  computed, 
leaving  a  6  x  2  matrix  rather  than  a  square  matrix.  The  term  —  is  defined  as 


Cl6c 

CLir 

CL6a 

Cy* € 

Cy6t 

Cy„ 

Cd6c 

Cd6t 

Cd6o 

Cl,. 

Cu, 

Ci6a 

Cmit 

Cmtr 

Cmta 

Cn(e 

cnt. 

Cn6a 
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Cfo  is  defined  to  be  the  vector  of  steady  state  coefficients: 


Cfo  = 


C  DO 

Cyo 

Clo 

ClO 
C mO 

Cn  0 


representing  conditions  in  trimmed,  balanced  flight.  This  definition  is  similar  to  the 
definition  used  by  Roskam  [Ref.  Ros  79].  In  other  references,  the  term  Cfo  can 
refer  to  the  nominal  value  of  the  coefficient  at  a  =  0.  However,  in  the  Taylor  series 
expansion  it  is  more  natural  to  use  the  first  definition  of  Cfo-'  therefore,  it  will  be 
used  in  the  following  derivation  and  modeling.  The  stability  and  control  derivatives 
are  usually  computed  in  the  so-called  wind  axis  coordinate  system.  The  wind  axis 
coordinate  system,  {IT},  is  defined  as  the  coordinate  system  that  results  when  the 
iB  axis  is  aligned  into  the  relative  wind.  This  axis  will  not  be  aligned  with  the  body 
coordinate  system  since  the  aircraft  flies  with  an  angle  of  attack,  a.  and  can  have  a 
sideslip  angle,  3.  The  transformation  from  {IT}  to  { B }  is  performed  in  the  same 
fashion  as  the  Euler  angle  transformations  mentioned  earlier.  The  rotation  matrix. 
B-R,  is  a  function  of  a  and  3.  and  is  expressed  as 


cos  q  cos  3  —  cos  q  sin  3  —  sin  o 
sin  3  cos  3  0 

sin  q  cos  3  —  sin  q  sin  J  cos  a 


(3.46) 


The  rotation  from  {IT}  to  { B }  is  made  by  premultiplying  the  force  or  moment  vector 
by  ®  R.  Additionally,  since  the  lift  and  drag  are  defined  as  positive  along  the  negative 
zb  and  /g  axes,  we  define  Faero  and  Naero  as 


(3.47) 


’  -D  1 

'  /  ‘ 

Faero  = 

Y 

-L 

and  A aero  = 

m 

n 
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In  order  to  write  Equation  3.45  in  state  space  form,  state  variables  must 
be  defined.  The  most  commonly  used  notation  to  use  for  the  state  vector  is  to  use 


x  = 


u 

v 

xv 

P 

<7 


(3.48) 


L  r  J 

However,  the  terms  x'  and  x'  in  Equations  3.40  and  3.41  cannot  be  used  directly  as 
states.  Define 


where 


and 


A/'  :  x  — >  x' 

M' :  x  —*  x' 

M'  =  diag{  1  /Vt,  1/Vj,  1/Vj,  6/21  r,  c/2VV,  6/2VV) 


(3.49) 


M'  = 


0  c/( 2Vt)  0  0  0  0 

0  0  6/(2 V 7)  0  0  0 


are  matrices  of  appropriate  dimensions.  The  complete  expression  for  B Faero  and 
B N aero  can  now  be  written  as 


B  Faero 


b 


a: 


AERO 


=  qS 


BR  0 


0 


B 

W 


R 


<Cf-°  +  f?w''  +  l>'  +  IA>  (3-50) 


and  can  be  substituted  into  Equation  3.38. 

2.  Other  Forces  and  Moments 

In  addition  to  the  forces  and  moments  due  the  aerodynamics  of  the  aircraft, 
forces  and  moments  due  to  the  propulsive  and  gravitational  forces  must  be  considered. 
Gravitational  forces  acting  on  the  aircraft,  B Fgrav,  can  be  found  by  premultiplying 
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1  Fgrav  by  the  appropriate  rotation  matrix,  where 


Fgrav  —  0 


Then 


»  T?  _  B  QL  r* 

*GRAV  —  irn  F GRAY  • 


(3.51) 


Propulsive  forces  and  moments,  B  Fprop  and  SA prop-  are  computed  directly  in  {B} 


and  can  be  expressed  as: 


'  Tx 

B  Fprop  =  Ty 

Tz 


» PROP 


= 


(3.53) 


where  T,'s  represent  the  forces  or  moments  due  to  thrust.  Computation  of  propulsive 
forces  and  moments  depends  on  each  particular  engine  installation,  and  must  be 
determined  for  the  individual  aircraft  modeled. 

Equations  3.51.  2.,  and  2.  can  now  be  substituted  into  Equation  3.3S: 


—  Bvbo 
dt  bujr 


b  ,  v 

—  WgX 


0  — S/fi1(S^’Bx  {bIbb^'b  +  IrB^'r))  J  [  b^'b 

f  ^  o  1  \  B1 


where 


0  Bln'  B\ 


bFgrav  bFprop 

0  B  N PROP 


\R  0 
0  IR 


?s{  Cro  +  gA/WgM'i  +  gA  }  |  }. 


In  order  to  write  Equation  3.38  in  state  space  form,  the  terms  associated  with  x'  must 


be  collected  and  moved  to  the  left  hand  side,  along  with  the  other  time  derivative 
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P  is  the  position  vector  of  the  aircraft,  and  A  is  the  matrix  of  kinematic  differential 
equations  based  on  Euler  angles  or  quaternions. 


IV.  COMPUTER  MODELING  OF  THE 
AIRCRAFT  EQUATIONS  O^  MOTION 


A.  IMPLEMENTATION  OF  EQUATIONS 
1.  Procedure 

The  AROD/Archytas  and  the  Bluebird  aircraft  models  require  different 
non-linear  equations  of  motion.  This  difference  is  due  to  the  unique  nature  of  each 
aircraft.  In  the  case  of  the  AROD.  the  angular  momentum  of  the  propeller  is  a  major 
factor  to  be  considered,  while  the  aerodynamics  effects  >n  a  hover  are  negligible. 
The  Bluebird  on  the  other  hand  is  a  conventional  aircraft,  and  exhibits  the  opposite 
characteristics.  Angular  momentum  from  the  propeller  is  small  and.  the  aircraft 
requires  the  stability  and  control  derivatives  associated  with  aerodynamic  flight.  All 
the  equations  were  implemented  in  a  systematic  manner  following  the  same  general 
approach  for  either  aircraft  model.  The  commercial  products  MatI.aB  and  SlMELINK. 
©  1990  1992  the  Mathworks.  were  chosen  for  the  modeling1,  mainly  due  to  the  ease 
of  expressing  matrix  equations.  The  model  was  constructed  using  the  following  steps. 

•  Kinematic  equations  of  motion  were  coded. 

•  Gravitational  forces,  with  the  direction  cosine  matrix  represented  by  Euler  an¬ 
gles  were  added  to  the  model 

•  Stability  and  control  derivatives  were  included  in  the  model,  as  well  as  engine 
thrust,  as  appropriate  for  the  aircraft  being  modeled. 


'All  rode  is  listed  in  APPENDIX  (’. 


•  Engine  and  actuators  were  added  to  the  model. 


•  Sensor  for  control  systems  were  added  to  the  model. 

In  the  first  stage  of  modeling,  analytic  linearization  was  also  carried  out  to 
verify  the  computer  calculations.  This  was  done  by  analytically  linearizing  the  matrix 
formed  from  the  six  non-linear  equations  governing  the  kinematic  motion  of  the  body. 
Nominal  values  wrere  substituted  in,  and  the  results  compared  to  the  trimmed  and 
linearized  values  obtained  from  the  SlMULINK  program2.  The  i.  xt  step  required  the 
addition  of  gravitational  terms,  and  analytic  linearization  was  still  manageable.  The 
linearized  matrix  now  included  nine  equations  and  nine  states  with  the  Euler  angle 
direction  cosine  matrix,  DCM,  and  ten  equations  with  ten  states  for  the  quaternion 
DCM.  Nominal  values  for  each  case  were  again  substituted  into  the  linearized  ma¬ 
trix  and  compared  to  the  linearized  model  derived  from  SlMULINK.  The  mclusion  of 
the  stability  and  control  derivatives  and  the  thrust  terms  presented  a  problem  that 
was  much  too  cumbersome  to  linearize  analytically.  Verification  of  the  data  at  this 
stage  was  accomplished  by  direct  comparison  of  the  dimensional  derivatives  resum¬ 
ing  from  numerical  linearization  of  the  plant.  In  the  case  of  Mie  AROD  the  results 
were  compared  with  data  published  by  Sandia  Labs  [Ref.  Wh  91];  for  the  Bluebird, 
eigenvalues  were  computed  and  then  compared  to  eigenvalues  obtained  by  classical 
analytic  methods  [Ref.  Sch  92]  and  [Ref.  Ros  79]. 

2.  AROD  Equations 

The  AROD  differs  from  the  Bluebird  primarily  in  that  the  aerodynamic 
forces  and  moments  are  negligible  while  the  craft  is  hovering.  The  powered  lift  does 
present  some  special  problems,  the  predominant  difficulty  being  the  gyroscopic  cou¬ 
pling  of  the  spinning  propeller.  Another  important  consideration  is  the  moment  due 

^Data  is  tabulated  in  APPENDIX  R. 
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to  the  swirl  effect  of  the  air  from  the  propeller  striking  the  vanes  mounted  aft  of  the 
propeller.  These  forces  and  moments  are  computed,  then  substituted  in  for  B F  and 
BN  in  Equation  3.38. 

a.  Control  Forces  and  Moments 

All  of  the  applied  forces  and  moments  in  the  AROD  are  due  to  the 
four  control  vanes  mounted  aft  of  the  propeller.  The  vanes  can  be  moved  in  different 
combinations,  as  discussed  in  later  in  this  chapter,  in  order  to  maneuver  the  AROD. 
The  forces  and  moments,  B F  and  B N ,  acting  on  the  AROD  are  computed  from  a 
Taylor  series  expansion  around  a  nominal  hover  point.  Since  aerodynamic  terms  are 
negligible, 

P Control 
P*  Control 

where 

.  a  =  [«„  e„  6.\ 

•  q,  =  1/2#>V;2 

•  V,  is  the  induced  velocity  through  the  propeller  [Ref.  Pro  90]  and 

V'2  =  T/(2Ap). 

•  A  is  the  inlet  area,  where  A  =  3.14  f2. 

•  R  is  the  propeller  radius,  where  R  =  1.0  f2. 

Notice  no  aerodynamic  forces  act  on  AROD  due  to  the  movement  of  the  elevator,  rud¬ 
der,  or  aileron  controls.  Again,  this  is  because  the  model  of  interest  is  only  designed 
to  fly  in  a  stable  hover.  The  other  forces  and  moments  involved  in  these  calculations 
are  due  to  gravitational  and  propulsive  action,  and  are  described  next. 


-  A  OdC  A 

q'ARd aa> 


(4.1) 
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b.  Additional  Forces  and  Moments 


The  gravitational  force  term,  bFgrav ,  is  computed  in  the  {B}  coor¬ 
dinate  system  with  the  rotation  matrix.  LBR.  For  the  AROD,  lBR  is  determined  by 
the  2-3-1  Euler  angle  rotation  sequence  defined  in  Equation  3.10  and  is  written  as 


B , 


Fgrav  =  BRL  Fgrav  where  L  Fgrav  — 


v 


0 
0 

L  m9 


(4.2) 


where 


B 


Fgrav  =  rn  g 


—  sin  0  cos 

sin©  sin  #  cos  <I>  +  sin  <*>  cos© 
-  sin  0  sin  'I'  sin  $  -f  cos  0  cos 


(4.3) 


The  forces  and  moments  due  to  the  propeller  thrust  were  determined 
experimentally  [Ref.  Sto  93]  and  are  discussed  in  a  later  section.  For  the  AROD,  the 
propulsive  force,  B Fprop-  is  acting  completely  along  the  xg  axis 


Fprop  = 


Ftx 

0 

0 


(4.4) 


where  Fjx  is  the  total  thrust,  determined  as  a  function  of  rpm. 

The  moment  resulting  from  thrust  is  SA prop  and  is  given  by  the 


vector  equation 


»N 


PROP  = 


W 

0 

0 


(4.5) 


where  lj  is  the  rolling  moment  due  to  thrust.  This  rolling  moment  is  due  to  the 
swirl  of  the  air  as  it  leaves  the  propeller  and  strikes  the  control  vanes.  This  term  was 
determined  experimentally  as  a  function  of  thrust,  and  is  also  discussed  further  in  a 
subsequent  section. 

When  all  the  terms  are  collected,  the  total  force  and  moment  applied 
to  AROD  can  be  expressed  as 


ejAR^AF 


B  Fprop 

B  A  PROP 


+ 


B 


Fgrav 

0 


(4.6) 
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The  force  and  moment  terms  in  Equation  4.6  can  be  substituted  into 
Equation  3.38  with  the  resulting  state  space  equation  given  by 


d_ 

Ti 

1  /m  0 
0 


B 


VBO 

3U>b 


HI- 


-bojb  x  bvBo 


rBx 


qAR^A  + 


B 1  y  {IbB^'B  +  IrP“JR) 

bFprop 
B  N prop 


+ 
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BF( 


GRA\ 
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J 

])) 


(4.7) 


Equation  4.7  can  now  be  written  in  the  state  space  form  and  programmed  using 
SlMULINK 


d 

Bi'bo 

_f 

—  bu>b  0 

dt 

BWJ3 

~l 

0  -fl/B-1(BWBx(B/B%  +  /flflWR) 

vBo 

?WB 
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I /m  0 
0 


1b 


B  Fprop 
BK> 


‘ PROP 
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B 


Fgrav 

0 


•  (4-8) 


3.  Model  Validation 

a.  Kinematic  Equations 

As  discussed  in  the  beginning  of  the  chapter,  the  first  step  in  the 
AROD  model  validation  was  to  linearize  the  non-linear  equations 


-j~Bvbo  =  1  /m(-BujB  x  Bi’bo) 


b 


dt 


WB 


Bj-\i  B.  ,  v  (t)  j  B  ,  .  B  j  B  ,  \ 

JB  ^9  x  1  Jr  u;r  +  Jr  wbJ, 


rB  j  _  B 


B 


(4.9) 

(4.10) 


with  the  result  given  by3: 


d_ 

dt 


6v 

<5w 

— w0  x  I'O  x 

Q  _B  J-\(B  r-B  ^  ,_s,\B  J_  I  B  , 


-bI-b1{b1rb^r  X  -(w0x)«/B  +  */SW0x) 


+ 


0 

—wq  x  B  Ir 


Sv 

6w 


SuJR.  (4.11) 


The  nominal  values  for  the  AROD  hover  operating  point, 

10 
0 
0 
0 
0 
0 


x  = 


(4.12) 


3See  APPENDIX  A  for  a  description  of  the  linearization  process 
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can  be  substituted  into  Equation  4.11  to  get 

'0000  0  O' 

0  0  0  0  0  10 

d_  6v  1  0  0  0  -10  0  0  \  bv  1 

dt  Su;  0  0  0  0  C  0  j 

0  0  0  0  0  —1.615 

.0  0  0  0  1.606  0 

where  the  linearized  matrix  associated  with  60,7*  is  zero,  since  —  x  B//?  =  0.  These 
values  match  the  values  obtained  by  linearizing  the  model  using  SlMl'LINK. 

The  inertial  cross  coupling  from  the  propeller  is  evident  from  the  lin¬ 
earized  results.  The  values  of  —1.6152  stc~l  and  1.6062  sec-1  are  defined  as  pitching 
moment  due  to  yaw  rate,  mT,  and  yawing  moment  due  to  pitch  rate.  nq.  respectively. 
The  gyroscopic  coupling  is  demonstrated  by  putting  a  step  moment  input  of  1  second 
duration  along  the  y  axis  and  observing  the  time  history  of  the  angular  rates  q  and 
r  as  shown  in  Figure  4.1.  The  AROD  has  a  tendency  to  spin  in  the  axis  orthogonal 
to  the  input  torque  as  shown  by  the  motion  r  about  the  z  axis. 


T«n«,  ttcondi 

Figure  4.1:  Gyroscopic  Motion  of  AROD 
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b.  Gravitational  Forets 


The  next  step  in  modeling  the  AROD  was  to  include  gravitational 
effects  into  the  model.  This  required  the  expression  for  B  Fc.ray  •  given  by  Equa¬ 
tion  4.3.  to  be  coded  into  a  MaTLAB  function.  This  function  block  was  then  added 
to  the  SlM U LINK  model,  as  shown  in  Figure  4.2.  The  equation  to  be  analyzed  at  this 


Figure  4.2:  Modeling  of  Gravitational  Effects  for  AROD 


stage  was 
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(4.14) 


where  5(A)  represents  the  kinematic  differential  equations,  defined  in  Equation  3.11. 
Notice  that  since  gravity  acts  through  the  c.g..  no  B Sgrav  term  exists.  Next.  Equa¬ 
tion  4.14  was  linearized  both  analytically  and  with  SlMl’LINK.  The  analytic  result. 


given  in  state  space  form,  was 


£ 

Jt 
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+  — *-^'0  x  B  Ir 

0 


(4.15) 


In  Equation  4.15,  the  functions  /(A),G(A),  and  h{ A.u,’)  are  partial  derivatives  of 
Equations  4.3  and  3.11  with  respect  to  A,  or  d/d\.  The  function  /(A)  is  the  linearized 


expression  for  B Fgrav  given  by  Equation  4.3,  where 


d  bFGrav 

dA  m 


/(A)  =  ^T-  —~=9- 


—  sin  0  cos  'P 
sin  0  sin  $  cos  4>  +  sin  $  cos  0 
—  sin  0  sin  $  sin  $  +  cos  0  cos  $ 


(4.16) 


—  COS0  cos'P 


sin0  sin'P 


g  —  sin0  sin'P  sin4>  +  cos4>  cos0  cos0  sin'P  cos$  —  sin$  sin0  sin0  cos'P  cos^ 

—  sin4>  cos0  —  sin©  sinvp  cos4>  —  cos4>  sin0  —  cos0  sin'P  sin4>  —  sin0  cos'P  sin4> 

(4.17) 

where  /(A)  is  evaluated  at  the  nominal  condition,  A0.  The  function 


G(A)  =  d/d*>(S(\)Bu;B 


is  derived  by  linearizing  the  kinematic  different •  ;1  equations  in  Equation  3.11.  The 

function  h(A.uj)  =  d/dA(S(A)B<DB)  is  not  presented  since  for  u,'0  =  0,  as  in  steady 

state  cruise,  h{ A,w)  =  0.  The  matrix  G'(A)  is  derived  from  the  following  equation 

q  q  1  —  cos'Ptan'P  sin4>tan^  p 

G{\)  =  —  (5(A)Bu.’b)  =  —  0  cos  $  sec  —  sin  sec  'P  q  (4.18) 

duJ  duJ  [  0  sin  $  cos  <J>  J  r 

1  —  cos  $  tan  $  sin  $  tan  $ 

0  cos  $  sec  —  sin  $  sec  $  .  (4.19) 

0  sin  $  cos  $ 
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These  matrices  were  assembled  as  in  Equation  4.15  and  were  evaluated  for  the  nominal 

flight  condition.  The  values  for  iq  for  this  flight  condition  were  given  by 

10  ' 

0 
0 
0 

0  .  (4.20) 

0 

0 
_7T 
2 

o  . 

After  substitution  into  Equation  4.15,  the  result  for  the  linearized  equations  of  motion 
was 

’0000  0  0  0  0  O' 

0000  0  10  0  0  32.174 

0  0  0  -10  0  0  0  -32.174  0 

,  [  6v  0000  0  0  0  0  0 

—  6-;  =  0  0  0  0  0  -1.615  0  0  0  (4.21) 

dt  <5A  J  0  0  0  0  1.606  0  0  0  0 

0001  0  00  0  0 

0000  1  00  0  0 

.0000  0  1  0  0  0. 

This  result  was  essentially  indentical  to  the  results  obtained  by  using  the  trim  and 
linearize  functions  in  SlMl  LINK.  The  gravitational  effects  acting  on  the  aircraft  were 
examined  by  running  a  simulation  of  the  non-linear  model  and  comparing  the  results 
to  those  from  Equation  4.11.  The  expectation  is  that  a  body  falling  in  earth's  gravity 
will  experience  an  acceleration  of  32.174//s2  and  as  shown  in  Figure  4.3,  this  expec¬ 
tation  is  realized  in  the  model.  In  Figure  4.3,  at  the  end  of  the  10  second  simulation, 
the  vertical  velocity,  u  in  this  case,  is  ~  320 //s,  as  was  predicted, 
c.  Additional  Forces  and  Moments 

The  last  step  in  modeling  the  AROD  was  to  include  the  forces  and 
moments  due  to  propulsive  and  control  action  that  act  on  the  aircraft.  The  data  used 
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0  2  4  6  8  10 

Time,  »*con<J» 

Figure  4.3:  Velocity  of  AROD,  Gravitational  Effects  Included 

in  modeling  these  control  forces  were  collected  experimentally  [Ref.  Sto  93]  and  then 
curve  fitted  to  a  function.  The  measurements  of  rolling  moment,  thrust,  rpm,  and 
vane  position  were  taken  for  various  configurations.  The  data  collected  were  then 
reduced  and  the  required  characteristics  were  computed.  The  accuracy  of  the  AROD 
model  depended  on  very  accurate  modeling  of  thrust  as  a  function  of  rpm,  moment  as 
a  function  of  thrust  produced,  and  moment  produced  by  deflecting  the  control  vanes 
in  the  different  combinations. 

Thrust  and  rolling  moments  were  measured  directly  at  different  power 
settings  ranging  from  3000  rpm  to  7600  rpm,  with  a  power  setting  of  ~  6400  rpm 
giving  a  thrust  of  ~  90/6/.  This  thrust  is  aprroximately  the  force  required  to  maintain 
a  hover  for  the  basic  AROD  configuration.  Figure  4.4  shows  the  linear  curve  fit 
through  the  thrust  vs.  rpm  data.  The  curve  was  fit  using  data  from  5000  rpm  to 
7600  rpm,  as  this  was  expected  to  approximate  the  normal  operation  range  in  flight. 
The  thrust  and  moment  data  are  plotted  in  Figure  4.5  along  with  the  line  fit  through 
those  data  points.  The  best  fit,  by  least  squares,  for  the  data  in  Figure  4.4  was  given 
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by 


Engtr*  Speed.  rpm 


Figure  4.4:  Thrust  vs.  RPM  for  AROD 


Figure  4.5:  Thrust  and  Moment  Data  for  AROD 


FTl  =  0.0297 Srvm  -  104.7,  (4.22) 


where  <5rpm  represents  the  rpm  at  a  given  throttle  setting.  The  best  fit  for  the  data 


in  Figure  4.5  was  given  by 


AV  =  —0.0542  Ft,  -  0.9138. 


(4.23) 


These  equations  were  then  used  in  modeling  the  engine  response  to  a  rpm  setting. 

The  engine  itself  was  modeled  using  a  second  order  transfer  function 
from  the  servo  position  to  6rpm ,  based  on  the  Sandia  Labs  [Ref.  Wh  87]  model.  The 
transfer  function  for  the  engine  was  given  as 


r  _  Ke 
rpm  S2  +  2  (u;nS  +  ^n 


6x 


(4.24) 


where  I\'e  =  900.  u>„  =  5 rad/ sec,  (  =  1.0,  and  6j  is  the  throttle  servo  position.  Since 
the  actual  throttle  position  is  set  via  a  radio  link  and  tests  have  not  been  set  up  to 
model  the  response,  it  was  ignored  in  the  model. 

The  engine  servo  could  be  modeled  and  a  transfer  function  from  com¬ 
manded  input  to  servo  position  was  determined.  This  transfer  function  also  was 
determined  by  Sandia  Labs  in  the  original  A  ROD  work  [Ref.  Wh  87]  to  be 

.2 


8j  — 


s'2  +  2(u.’ns  +  u)2 


Uj 


(4.25) 


where  uin  =  20 rad./ sec,  (  =  0.6,  and  uj  represents  the  commanded  input  to  the  servo. 
It  should  be  noted  that  Equation  4.25  was  also  used  for  the  control  vane  servos,  since 
all  the  servos  in  the  AROD  were  identical. 

In  order  to  model  the  moments  due  to  control  surface  deflection, 
B ^control,  computation  of  control  vane  effectiveness  was  necessary.  Again,  the  data 
gathered  by  [Ref.  Sto  93]  was  used  to  compute  vane  effectiveness  for  the  AROD 
model.  Figure  4.6  shows  the  moment  data  for  75%  thrust.  This  power  setting  cor¬ 
responds  closely  to  the  thrust  required  to  maintain  hovering  flight.  A  line  was  fitted 
through  the  data  points,  and  the  slope  of  this  line  was  used  to  determine  the  rolling 
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Figure  4.6:  Moment  Due  to  Vane  Deflection 


moment,/,  for  a  given  vane  deflection.  The  equation  for  the  curve  was 

/  =  0.5523  VAvg  ~  5.7248  (4.26) 


for  6580  rpm,  where  Vavg  was  the  average  vane  deflection,  1/4  \  j.  in  degrees. 

This  averaging  was  required  since  the  individual  vanes  were  at  slightly  different  posi¬ 
tions  with  respect  to  the  commanded  position.  Computing  the  vane  effectiveness  for 
roll  was  done  by  differentiating  with  respect  to  the  vane  deflection. 


h.  = 


dl 


=  0.5523 


ft  -  Ibj 


(4.27 


OVavg  ~  d(9 

The  rolling  moment  should  be  non-dimensionalized  for  use  in  the  equations  of  motion 
given  by  Equation  4.8.  This  was  done  by  applying  the  following  equation 


Cn  = 


h 


1/2  p\'2Sb' 


(4.28) 


To  use  the  measurable  quantities  available,  redefining  the  terms  in  Equation  4.28  was 
necessary.  The  representative  area,  S.  is  defined  as  the  inlet  area  to  the  duct.  A.  The 
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characteristic  length,  b,  is  defined  as  the  propeller  radius.  R.  The  velocity  term  was 
defined  as  induced  velocity,  \ .  With  these  terms  the  non-dimensional  coefficient  C/# 
could  be  calculated  as 

c".  -  Tmtm  (4-29> 

and  was  computed  to  be  Ci6a  =  1.438/rad.  The  rolling  moment  was  measured  during 
the  testing  and  was  easily  non-dimensionalized  in  a  form  suited  to  computer  modeling. 
However,  no  measured  data  existed  for  the  pitching  and  yawing  moments.  This 
required  the  estimation  of  pitch  and  yaw  vane  effectiveness  by  a  simple  ratio  technique. 
First,  it  was  assumed  that  the  vane  effectiveness.  lsa  for  two  vanes  (2V7),  would  be 
exactly  twice  the  /$a  for  one  vane,  and  1/2  the  lsa  for  four  vanes,  (4V).  This  would 
make 


l-Jt—hv  =  hJF-Uv  (4.30) 


dV.y 


AVG 


2'3V 


AVG 


or  numerically, 


(JU,  =  0  W'-“' 


'dVAva”'  deg 

Now,  the  ratio  of  the  dimensional  derivatives  to  the  moment  arms  was  set  up  as 

(  91  \  /  dm 

[  9VAVG  >2V  _  l  avAV 

L  l 


(4.31) 


:) 


(4.32) 


where  1T  was  the  distance  from  the  c.g.  along  the  r-axis.  to  the  midspan  of  the  control 
vane,  a  distance  of  9.0  in.  The  distance  ly  was  the  distance  from  the  c.g.  along  the 
y-axis,  to  the  1/4  chord  of  the  control  vane,  a  distance  of  15.43  in.  The  calculation 
was  performed  and  the  result  for  the  pitching  moment  was 


dm  ft  —  Ibj 

=  0.4735- - 1 


(4.33) 


dV avg  deg 

This  was  non-dimensionalized  using  Equation  4.28  with  the  result  of  Cm(e  =  1.2332 rad-1. 
Moreover  because  of  the  symmetrical  design  of  AROD,  Cm(f  =  Cnt  .  The  results  for 


41 


vane  effectiveness  are  important  to  a  high  fidelity  model  and  are  presented  along 
with  other  relevant  data  in  Table  4.1.  The  non-dimensional  derivatives  in  Table  4.1 

TABLE  4.1:  NON-DIMENSIONAL  DERIVATIVE  DATA  FOR  AROD 


IMB 

rad  1 

deg~l 

Positive  Vane  Deflection 

113 

1.438 

0.0251 

V'l  +  V'2  +  V3  +  1  4 

mm 

-1.233 

-0.0215 

r4  -  v  2 

IBM 

-1.233 

-0.0215 

is  -  r, 

were  substituted  into  the  term  ( dCF/dA )  in  Equation  4.7.  Written  as  a  matrix,  the 
non-dimensional  derivatives  are 


dCF 


A  = 


'cmtt  0  0  ' 

'  St  ' 

0  0 

br 

1 

0 

0 

e 

1 _ 

K 

(4.34) 


Equations  4.22,  4.23,  and  4.34  were  added  as  a  separate  block  in  the  model,  resulting 
in  a  block  diagram  shown  in  Figure  4.7. 


Figure  4.7:  Non-Linear  Equations  of  Motion  Model 


In  the  SlMl'LINK  linearization  results,  no  change  was  expected  in  the 


A  matrix  since  no  aerodynamic  terms  were  added  to  the  model.  Rolling  motion  is 
shown  in  Figure  4.8  as  an  example  of  the  unstable  natu~e  of  the  AROD  without  a 
stability  augmentation  system. 


Figure  4.8:  Rolling  Motion  For  Complete  Non-Linear  Equations 


4.  Bluebird  Equations  of  Motion 

a.  hint  mafic  Equations 

Modeling  the  kinematic  equations  of  motion  for  Bluebird  was  accom¬ 
plished  using  the  same  procedure  as  was  used  in  the  AROD  modeling.  A  slight  dif¬ 
ference  arose  in  the  derivation  of  the  angular  acceleration  equation  for  Bluebird  since 
the  term,  Irujb,  for  the  angular  momentum  due  to  propeller  rotation  is  negligible. 
Thus,  the  following  equations  were  used  in  the  first  stage  of  modeling. 

-j- Bvb°  —  ~B^'b  x  Bvbo  (4.35) 

~77 B x  BIBBojB).  (4.3b) 
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These  equations  were  coded  into  a  MaTLAB  function  and  placed  into  a  block  diagram, 
shown  in  Figure  -1  9.  Equations  4.35  and  4.36  w-ere  linearized  analytically'1  with  the 


lanibdadol 


Figure  4.9:  Block  Diagram  of  Kinematic  Equations  of  Motion 


resulting  state  space  equation 

d_  6v  _  u>0  x  t’ox  6v 

d/[6u;J_[  0  -B/e1{-(u;ox)B/fi  +  B/B(^oX)} 

This  equation  was  evaluated  at  the  nominal  flight  conditions,  determined  by  the 
typical  cruise  condition  for  the  Bluebird  aircraft.  A  state  vector  for  the  nominal 
flight  conditions  is  given  by 

72.9954  f/s 

0 

6.6757  f/s 
0 
0 
0 

“See  APPENDIX  A  for  the  details 


(4.38) 
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These  values  of  Jo  were  substituted  into  Equation  4.37: 


'0  0  0  0  -6.675  0 

0  0  0  6.675  0  -72.995 

d_  bv  _  0  0  0  0  72.995  0 

dt  6*j  _  0  0  0  0  0  0 

0  0  0  0  0  0 

.0  0  0  0  0  0 

These  results  were  in  complete  agreement  with  the  data  obtained  by  trimming  and 
linearizing  the  non-linear  model  with  the  SlMl  LINK  program.  The  comparison  be¬ 
tween  Equation  4.39  and  Equation  4.13  shows  the  absence  of  any  angular  rate  cross 
coupling  in  Bluebird.  The  absence  of  the  cross  coupling  terms  results  from  the  choice 
to  ignore  the  angular  momentum  contribution  from  the  propeller, 
b.  Gravitational  Forces 

After  the  basic  kinematic  equations  Equation  4.35  and  4.36  were  put 
into  block  diagram  form,  it  was  an  easy  matter  to  include  additional  blocks.  The 
next  block  was  a  function  block  including  the  B Fgray  terms.  The  model  with  the 
gravitational  terms  included  is  shown  in  Figure  4.10.  The  equations  of  motion  to  be 
modeled  at  this  stage  were 


(4.40) 

(4.41) 

(4.42) 


where  5(A)  is  the  set  of  kinematic  differential  equations  based  on  the  3-2-1  Euler 
angle  rotation  in  Equation  3.7.  These  equations  were  then  linearized  analytically, 
with  the  result 


,  [  6v  ’  -u;0x  t’ox  (//(A)  I  ’  8v 

y  6*'  =  0  -s/e1{-(-’ox)s/e-f  S/B(^0x)}  0  ^ 

dt  [  6\  \  [  0  (7(A)  A(A,-:)  J  5\  . 
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(4.43) 


Figure  4.10:  Gravitational  Forces  Model 


where  the  terms  /(A),  G(A),  and  h( A,u>)  are  based  on  the  3-2-1  Euler  rotation  se¬ 
quence.  The  linearization  is  done  in  the  same  manner  as  was  shown  in  Equations  4.16. 
4.17,  and  4.19  resulting  in 

0  —  cos  0  0 

/(A)  =  cos0cos$  —  sin0sin$  0  ,  (4.44) 

—  cos  0  cos  $  —  sin  0  cos  $  0 


and 


Note  that 


1  sin4>tan0  cos  $  tan  0 
G( A)  =  0  cos$  —  sin  4> 

0  sin  $  sec  0  cos<I>sec0 

k( A,  w)  =  ^5(A)b^|0  =  0, 


(4.45) 
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since  u?o  =  0.  When  Equation  4.43  is  evaluated  at  the  nominal  condition  Jo  given  by 

’  72.9954  f/s  ' 

0 

6.6757  f/s 
0 

xo  =  0 

0 
0 

0.0912  rad 
0 

the  resulting  equation  is 

'  0  0  0  0  -6.675  0  0  -32.043  0  ' 

0  0  0  6.675  0  -72.995  32.043  0  0 

0  0  0  0  72.995  0  0  -2.930  0 

,  6v  1  0  0  0  0  0  0  0  0  0  6v  ' 

—  6u;  =  0  0  0  0  0  0  0  0  0  6u  . 

0000  0  0  0  0  0  [  6A 

0  0  0  1  0  0.0915  0  0  0 

0000  1  0  0  00 

.0  0  0  0  0  1.004  0  0  0. 

(4.46) 

In  this  instance,  as  well,  the  results  of  the  analytic  linearization  in  Equation  4.46 
agree  very  closely  with  the  computed  results. 

A  non-linear  simulation  of  the  system  in  Figure  4.10  should  show  an 
increasing  downward  velocity  due  to  the  gravitational  effects  of  B  Fgrav-  One  would 
also  expect  to  see  decreasing  forward  velocity  due  to  the  “drag-like’'  term  that  arises 
with  the  introduction  of  the  angle,  O.  This  plot  is  shown  in  Figure  4.11. 
c.  Aerodynamic  Forces  and  Moments 

Completion  of  the  Bluebird  equations  of  motion  model  required  the 
modeling  of  the  aerodynamic  forces  and  moments  acting  on  the  aircraft.  A  simple 
engine  model  was  also  developed.  No  analytic  linearization  was  performed  at  this 
stage  due  to  the  increased  complexity  of  the  model.  Verification  of  the  computer 
results  was  accomplished  by  comparing  the  modes  and  eigenvalues  of  the  computed 


Figure  4.11:  Gravitational  Effects  on  Velocity 


plant  with  those  resulting  from  substituting  the  stability  and  control  derivatives  into 
equations  developed  in  [Ref.  Sch  92]. 

The  aerodynamic  forces  and  moments  as  described  in  Equation  3.50 
were  coded  as  a  MaTLAB  function,  then  included  as  a  block  in  the  model  shown  in 
Figure  4.12. 

Next,  it  was  necessary  to  premultiply  all  the  blocks  by  \_1,  since 
this  added  the  effects  of  the  a  derivatives.  Now  the  important  task  was  to  calculate 
the  stability  and  control  derivatives  using  the  general  reference  for  the  estimation  of 
non-dimensional  derivatives,  DATCOM  [Ref.  USAF  60].  The  stability  and  control 
derivatives  wrere  computed  based  the  aircraft  geometry  and  the  control  surfaces.  These 
values  are  tabulated  in  Table  4.2  and  in  Table  4.3  where  the  non-dimensional  force  is 
listed  on  the  left  side,  and  the  particular  derivative  is  determined  using  the  top  row. 
For  example,  Coa  =0.188  using  Table  4.2.  The  terms  in  CV0  were  also  estimated  to 
be  Cl0  =  0.385  and  Cp0  ~  0.03,  with  all  other  terms  equal  to  zero  since  the  aircraft 
is  in  straight  and  level  flight  at  this  trim  condition. 
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Figure  4.12:  Full  Non-Linear  Equations  of  Motion  Model 
TABLE  4.2:  NON-DIMENSIONAL  STABILITY  DERIVATIVES 


D 

a 

P 

q 

r 

d 

Cd 

Q 

0 

0.188 

0 

0 

0 

0 

Cy 

D 

-0.31 

0 

0 

0 

0.0973 

0 

Cl 

D 

0 

4.22 

0 

3.94 

0 

1.32 

G 

Q 

-0.0597 

0 

-0.363 

0 

0.100 

0 

cm 

Q 

0 

-1.163 

0 

-11.77 

0 

-4.70 

IBM 

a 

0.0487 

0 

-0.0481 

0 

-0.0452 

0 

The  B Fpfiop  was  based  on  estimating  engine  thrust  for  a  4  HP  engine 
and  a  propeller  efficency,  pp ,  of  0.65.  Thrust,  To  was  estimated  using  the  equation 


r=  550 n^HPp,  (4.47) 

b'O  Po 

where  pn  is  the  density  at  the  operating  altitude,  p0  is  the  sea  level  density,  and  Co  is 
the  velocity  of  the  aircraft  in  i/s.  The  result  was  an  engine  thrust  of  T0  =  19.5/6/,  for 
a  density  ratio  of  1.  This  could  be  directly  factored  into  the  equations  of  motion  as 
7o<5j,  where  <5j  was  arbitrarily  scaled  from  zero  to  one.  8j  =  0  represents  zero  thrust 
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TABLE  4.3:  NON-DIMENSIONAL  CONTROL  DERIVATIVES 


8e 

Sr 

Sa 

Co 

0.065 

0 

0 

Cy 

0 

0.0697 

0 

CL 

0.472 

0 

0 

Ci 

0 

0.0028 

0.265 

cm 

-1.410 

0 

0 

cn 

0 

-0.0329 

-0.0347 

and  6j  =  1  represents  maximum  thrust. 

The  preceeding  values  were  substituted  into  the  appropriate  MaTLAB 
functions  and  the  entire  equations  of  motion  model  was  trimmed,  linearized,  and 
then  compared  with  the  analytic  results  based  on  classical  techniques  for  determining 
eigenvalues  and  eigenvectors.  The  eigenvalues  are  compared  in  Table  4.4.  The  results 
from  the  computed  eigenvalues  is  very  close  to  the  eigenvalues  derived  by  analytic 
methods. 

TABLE  4.4:  COMPARISON  OF  EIGENVALUES  FOR  BLUEBIRD 


MODE 

COMPUTED 

ANALYTIC 

LONGITUDINAL 

Phugoid 

-0.0191  ±0.4963; 

-0.0473  +  0.4940; 

Short  Period 

-3.9833  ±  3.5521 ;' 

-4.0034  ±  3.5462; 

LATERAL 

Dutch  Roll 

-0.522  +  3.6194; 

Short  Period 

-5.5629 

-5.6654 

Spiral 

+0.0420 

+0.0420 

B.  VALIDATION  OF  AN  INDEPENDENT  CASE 

Although  verification  of  the  model  was  accomplished  at  each  stage,  comparison 
of  the  results  from  the  numerical  linearization  with  linearized  results  from  an  inde- 
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pendent  source  was  still  necessary  before  the  model  could  be  considered  completely 
reliable.  The  test  case  selected  was  the  Cessna  172  documented  in  [Ref.  Ros  79]  as 
airplane  A.  The  tabulated  non-dimensional  derivatives  and  given  flight  conditions 
were  used  as  inputs  to  the  non-linear  model.  The  model  was  then  trimmed  at  the 
specified  flight  condition  of  Vj  =  219  f  js  and  0o  =  0.  Using  the  resulting  state  and 
input  vector,  the  model  was  linearized  around  thesed  nominal  conditions.  With  the 
linearized  plant,  eigenvalues  were  determined  and  compared  to  those  tabulated  for 
airplane  A  (see  Table  4.5).  Very  little  difference  between  the  eigenvalues  is  seen  in 
the  longitudinal  modes.  Slightly  greater  differences  are  noticed  when  comparing  the 
lateral  modes,  but  these  differences  are  still  fairly  small.  Another  very  good  method 

TABLE  4.5:  EIGENVALUE  COMPARISON  FOR  CESSNA  172  TEST 
CASE 


Mode 

Numerical 

Independent 

Longitudinal 

Short  Period 

-4.130  ±4.3895; 

-4.130  ±4.390; 

Phugoid 

-0.0209  ±  0.1794; 

-0.02092  ±  0.1797; 

Lateral-Directional 

Dutch  Roll 

-0.6947  ±  3.3080; 

-0.6858  ±3.306; 

Spiral 

-12.4309 

-12.43 

Roll  Response 

-0.0109 

-0.01095 

for  comparing  the  results  of  the  numerical  linearization  with  RoskanTs  tabulated  data 
is  to  form  plant  and  control  matrices  for  the  test  aircraft,  using  the  linear  algebraic 
method  taught  in  AE  3340  [Ref.  Sch  92].  The  resulting  A  and  B  matrices  were  then 
compared  by  applying  step  elevator,  rudder,  and  aileron  inputs  and  plotting  the  re¬ 
sults.  The  results  for  the  step  elevator  input  are  given  in  Figure  4.13.  These  plots 
show  very  little  difference  in  the  vertical  velocity,  w.  Differing  amplitudes  are  shown 
for  the  horizontal  velocity,  but  since  the  non-dimensional  velocity,  u / Vj ,  from  the 
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dimensional  derivatives  was  scaled  to  be  equivalent  to  the  state,  u.  computed  in  the 
numerical  linearization,  these  magnitude  errors  are  not  indicative  of  a  poor  model, 
only  a  slight  difference  in  the  computed  damping  is  shown.  The  results  from  the 
analytic  model  were  scaled  up  by  the  nominal  airspeed,  Vj  to  compare  with  the  nu¬ 
merically  linearized  results.  This  would  have  the  effect  of  magnifying  any  differences 
between  the  analytic  model  and  the  numerical  model.  The  natural  frequency  of  both 
the  analytic  and  numerical  results  are  quite  similar. 


Figure  4.13:  Comparison  of  Longitudinal  Responses  to  Step  Elevator  Input 

Lateral  responses  to  step  rudder  and  step  aileron  inputs  are  shown  in  Fig¬ 
ures  4.14  and  4.15.  Very  little  difference  between  the  models  is  visible  in  these  plots. 
The  errors  are  shown  in  Figures  4.16  -  4.1& 

It  can  be  concluded  that  the  results  from  the  numerical  linearization  are  quite 
close  and  furthermore  that  the  equations  used  in  the  model  are  indeed  correct.  More¬ 
over,  the  linearized  results  presented  for  both  the  A  ROD  and  Bluebird  aircraft  should 
be  considered  accurate  and  are  suitable  for  the  Guidance,  Control,  and  Navigation 
system  designs  that  will  follow. 
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Plot  of  Stop  Rudder  Dots 


Time  teconds 

Figure  4.14:  Comparison  of  Lateral  Responses  to  Step  Rudder  Input 


Tim*.  Mconds 


Figure  4.15:  Comparison  of  Lateral  Responses  to  Step  Aileron  Input 
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Figure  4.16:  Difference  in  Analytic  and  Numerical  Results,  Step  Elevator 
Input 


Figure  4.17:  Difference  in  Analytic  and  Numerical  Results,  Step  Rudder 
Input 
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Figure  4.18:  Difference  in  Analytic  and  Numerical  Results, 
Input 
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V.  SENSOR  AND  ACTUATOR  MODELING 

The  AROD  engine  and  actuators  were  modeled  as  a  second  order  transfer  func¬ 
tions,  based  on  data  collected  by  Sandia  Labs  [Ref.  Wh  87].  Sensors  were  also  mod¬ 
eled  as  second  order  transfer  function  based  on  dat  ^  supplied  by  Watson  Industries 
[Ref.  WAT  93]. 

The  complete  non-linear  model  for  an  aircraft  shou'd  include  models  of  the 
sensors  on  board  for  measurement  of  acceleration,  angular  rates,  pitch  and  roll  an¬ 
gles,  and  headings.  The  inertial  device  chosen  for  the  AROD  project  was  the  Watson 
Industries  IMU-600D  inertial  measuring  unit.  This  device  contains  a  triaxial  ac¬ 
celerometer.  a  triaxial  rate  sensor,  two  liquid  pendulous  devices  (for  bank  and  pitch 
angle),  and  a  magnetic  heading  indicator.  The  characteristics  of  these  devices  must 
be  ac'  mtely  modeled,  since  it  is  the  sensor  output  that  drives  the  control  system. 
In  the  rest  of  this  section,  the  accelerometers,  rate  gyros,  and  inclinometers  will  be 
modeled,  as  well  as  the  sources  of  error  inherent  in  the  design  of  the  sensor  devices. 

A.  ACCELEROMETER  MODELING 

The  term  accelerometer  is  not  entirely  accurate,  since  the  device  does  not  mea¬ 
sure  true  acceleration,  but  rather  the  difference  between  acceleration  and  gravity 
[Ref.  Bro  64].  This  effect  is  referred  to  as  the  Einstein  Uncertainty  Principle ,  and  is 
represented  in  equation  form  as 

f  =  g-a  (5.1) 

where  g  and  a  are  the  specific  forces  of  gravity  and  acceleration  of  the  aircraft  and 
are  measured  positive  downwards.  The  accelerometer  mode!  relevant  to  this  equation 
is  shown  in  Figure  5.1.  The  tri  axial  accelerometer  of  the  IMT-60UD  can  be  modeled 
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Figure  5.1:  Typical  Accelerometer  Model 


as  three  simple  single-axis  accelerometers,  as  has  been  established  through  conversa¬ 
tions  with  the  manufacturer  [Ref.  WAT  93].  A  schematic  representation  basic  device 
pictured  in  Figure  5.1  is  modeled  by  an  ordinary  differential  equation  [Ref.  Sil  91] 


..  0  .  k 

x  H - x  +  — x  =  -y 

m  m 


(5.2) 


where  x  denotes  the  displacement  of  the  mass  from  its  equilibrium  position  and 
y  —  g  —  ci  is  the  projection  along  the  case  axis  of  the  vector  sum  of  gravity  and 
acceleration.  The  terms.  0,  k,  and  m  represent  the  damping,  spring  coefficient,  and 
mass,  respectively,  of  the  device.  The  accelerometer  described  in  Equation  o. 2  can 
be  modeled  as  a  second  order  low  pass  filter,  but  the  actual  accelerometer  has  a  flat 
response  up  to  1000Hz,  so  it  was  not  modeled.  A  thirl  order  Chebychev  anti-aliasing 
filter  with  a  cut-off  frequency  of  20  Hz  was  added  to  the  accelerometers.  This  filter 
is  the  device  that  was  modeled.  The  Chebychev  filter  gives  the  advantage  of  a  flat 
passband,  and  a  very  sharp  drop  off  at  the  cut-off  frequency.  The  equation  used  to 
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describe  a  Chebychev  filter  [Ref.  St  88]  is 


Wlp(^)\2 


1 


(5.3) 


where  C/v  is  the  Nth  order  Chebychev  polynomial,  t  is  the  parameter  thats  sets  the 
ripple  in  the  passband,  and  \Hpp(ju;)\2  is  the  magnitude  of  the  filter.  It  was  not 
necessary  to  compute  the  filter  analytically,  as  SlMULINK  provides  a  block  function 
which  performs  the  required  steps  based  on  the  passband  ripple  of  0.1  db  and  the 
desired  cutoff  frequency  of  20  Hz.  The  block  diagram  for  the  accelerometer  model 
is  shown  in  Figure  5.2.  Included  in  this  diagram  are  the  blocks  containing  the  error 
modeling,  as  well  as  the  blocks  used  to  correct  for  the  Einstein  uncertainty.  Figure  5.3 
shows  the  linear  synthesis  model  used  to  generate  the  accelerations  to  drive  the  sensor 
models.  The  synthesis  model  was  derived  from  the  Bluebird  model  discussed  in 
Chapter  IV. 


Figure  5.2:  Accelerometer  Modeling 
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Figure  5.3:  Synthesis  Model  for  Accelerometers 
1.  Error  Model 

No  matter  how  well  the  sensor  device  is  constructed,  any  accelerometer 
is  subject  to  certain  errors  in  the  linear  acceleration  measurements.  These  errors 
can  occur  for  several  reasons;  some  as  mechanical  and  others  are  due  to  the  physical 
placement  of  the  accelerometer  on  the  aircraft.  The  mechanical  errors  accounted  for 
in  the  IMU-600D  tri-axial  accelerometers  are 


•  Acceleration  Bias 

-  average  readings  for  -fig  and  -1  g  loads 


•  Acceleration  Scale  Factor  error 

-average  difference  between  readings  from  -f  1  g  and  -\g  loads 
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•  Cross  Axis  Sensitivity  errors 

-measurements  due  to  the  misalignment  of  an  accelerometer  with  the  appropri¬ 
ate  axis 

-measurements  of  off-axis  accelerations  are  measured 

•  Acceleration  Noise  Floor 

-threshold  below'  which  measurements  of  acceleration  can  not  be  made 

Errors  in  the  measured  accelerations  can  also  occur  if  the  accelerometers 
are  not  located  at  the  c.g.  of  the  aircraft,  since  arbitrarily  located  points  on  a  body 
will  experience  additional  accelerations  due  the  the  angular  momentum  of  the  body. 
A  mathematical  correction  for  the  error  will  be  presented  later,  and  will  describe  the 
angular  and  linear  accelerations  of  an  arbitrary  point  on  a  rigid  body.  However,  the 
correction  will  not  be  applied  to  the  models  here,  since  the  correction  is  expected  to 
have  a  very  small  effect  on  the  sensor  measurements  due  to  the  small  displacements 
away  from  the  aircraft  c.g. 

Error  terms  are  quantified  in  terms  of  full-scale  measurements.  The  me¬ 
chanical  errors  are  tabulated  in  Table  5.1  along  with  other  important  characteristics. 
Figure  5.4  shows  the  block  diagram  with  error  inputs  applied  to  the  accelerations 

TABLE  5.1:  ACCELEROMETER  CHARACTERISTICS 


Acceleration  Range 

±2g's 

Acceleration  Bandwidth 

20  Hz 

Acceleration  Bias 

0.2%  of  Full  Scale 

Acceleration  Scale  Factor 

0.2%  of  Full  Scale 

Acceleration  Noise  Floor 

0.0005  g's 

Cross  Axis  Sensitivity 

0.5%  of  Full  Scale 
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Figure  5.4:  Error  Model  for  Accelerometers 

measured  by  the  accelerometers.  The  output  from  the  error  computations  in  Fig¬ 
ure  5.4  is  the  measured  acceleration  from  the  accelerometer  output  to  the  control 
system.  The  cross  axis  terms  are  determined  through  a  matrix 


x  = 


0  ty  t: 

1 1  0  e. 

tX  ty  0 


X 


(5.4) 


where  e,  is  the  cross  axis  error  term  and  x  is  the  error  in  the  acceleration  due  to  the 
cross  axis  sensitivity. 

2.  Results  and  Validation 

Accelerations  were  measured  for  step  aileron,  elevator,  and  rudder  inputs, 
and  the  results  then  compared  to  the  actual  accelerations  computed  for  the  equations 
of  motion  for  the  aircraft.  A  linear  synthesis  model  was  used  for  the  initial  testing. 
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while  the  results  from  the  non-linear  mode!  were  used  for  validation  of  the  accelerom¬ 
eter  model.  Figures  5.5,  5.6,  and  5.7  show  comparisons  between  measured  and  actual 
accelerations  generated  from  simulations  of  the  non-linear  model.  The  accelerations 
computed  for  the  longitudinal  and  lateral  cases  were  in  close  agreement  with  the 
accelerations  computed  by  the  non-linear  model. 


Figure  5.5:  Measured  and  True  Acceleration  From  a  Step  Elevator  Input 
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Acceleration,  g's  T)  Acceleration,  g’s 


igure  5.6:  Measured  and  True  Acceleration  From  a  Step  Aileron  Input 


Figure  5.7:  Measured  and  True  Acceleration  From  a  Step  Rudder  Input 
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B.  RATE  GYRO  MODELING 


The  rate  gyros  consist  of  a  rotating  disk,  mounted  in  a  gimbal  mechanism. 
Though  both  single  and  two  degree  of  freedom  gyros  are  common,  the  tri-axial  angle 
rate  gyro  supplied  in  the  IMU-600D  is  modeled  as  three  single-dcgree-of-freedom 
rate  gyros,  each  measuring  the  angular  rate  along  a  particular  axis.  The  dynamics  of 
a  gyro  can  be  modeled  [Ref.  Sil  91]  using  Euler's  law 


B 


No  = 


B  L 


G 


wher  {B}  and  {6’}  are  the  body  and  gyro  coordinate  systems,  respectively.  This  can 
be  expanded  to 


V  _  B  ,  B  /  ,  B  ua  G  j 

■N°  -  Lc  +  “HI,  La 


where  G L a  =  IaB^'a  and  the  time  derivative.  j^Lg  is  zero  when  the  wheel  rotates 
at  a  constant  speed.  Thus,  except  for  the  period  when  the  wheel  is  coming  up  to 
speed,  the  equation  for  gyroscopic  motion  in  a  rate  gyro  can  be  written  as 


AV;  = 


•G 


(5.5) 


where  BXa  is  the  torque  vector  acting  on  the  gyro  element  and  B^a  is  the  angular 
rate  of  the  gyro  frame.  Transfer  functions  can  be  developed  for  the  torque  input  to 
the  input  axis  as  shown  in  Figure  5.8  and  the  rate  output  of  the  output  axis.  This 
derivation  was  not  developed  since  data  required  for  the  gyro  disk  inert -a,  la  and  the 
speed  of  rotation,  bojg,  among  other  terms,  is  not  available  from  the  manufacturer. 
The  rate  gyros  were  modeled  as  second  order  transfer  functions,  with  u^n  =  50  Hz. 
A  third  order  Chebychev  filter  with  a  cut-off  frequency  of  20  Hz  was  added  to  the 
gyro  model  as  an  anti-aliasing  filter.  This  Chebyshev  filter  was  identical  to  the  ones 
described  for  the  accelerometers.  The  block  diagram  of  the  rate  gyro  model  is  shown 
in  Figure  5.9.  The  linear  synthesis  model  is  shown  in  Figure  5.10  and  was  used  to 
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test  the  gyro  models.  In  Figure  5.9  the  gyro  elements  are  shown  for  each  axis  and 
blocks  for  the  error  calculations  are  shown  as  well. 

1.  Error  Modeling 

Error  terms  are  also  present  in  the  rate  gyros  and  are  due  to  either  physical 
location  on  the  aircraft,  or  mechanical  errors,  as  in  the  accelerometers.  Errors  due 
to  physical  placement  away  from  the  c.g.  can  be  corrected  by  using  the  equations 
derived  for  the  linear  and  angular  acceleration  of  an  arbitrary  point  on  the  aircraft, 
shown  later  in  this  chapter.  Mechanical  errors  are  defined  in  the  same  manner  as  was 
done  for  accelerometers  in  Section  A,  and  are  listed  in  Table  5.2  along  with  other 
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TABLE  5.2:  GYRO  CHARACTERISTICS 


Rotational  Rate  Range 

±114.6de<7/sec 

Rotational  Rate  Bandwidth 

20  Hz 

Rotational  Rate  Bias 

2.0%  of  Full  Scale 

Rotational  Rate  Scale  Factor 

0.5%  of  Full  Scale 

Rotational  Rate  Noise  Floor 

0.05%  of  Full  Scale 

Cross  Axis  Sensitivity 

0.5%)  of  Full  Scale 

important  characteristics.  The  cross  axis  error  is  modeled  in  the  same  way  as  for 
the  accelerometers  by  the  use  of  the  same  matrix  for  computing  errors  in  the  angular 
rates. 

2.  Results  and  Validation 

Angular  rates  were  measured  for  step  aileron,  elevator,  and  rudder  inputs 
and  compared  to  the  actual  angular  rates  computed  from  the  equations  of  motion. 
First  a  linearized  synthesis  model  was  used  in  the  testing  stages;  then  the  sensors 
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Figure  5.10:  Rate  Gyro  Synthesis  Model 

were  integrated  with  the  non-linear  aircraft  simulation.  Comparisons  of  actual  and 
measured  acceleration  are  given  in  Figures  5.11,  5.12,  and  5.13.  These  figures  show 
the  accuracy  of  the  angular  rate  sensors  in  measuring  the  angular  rates  computed 
with  the  non-linear  equation  of  motion  model. 
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Figure  5.11:  Measured  and  True  Angular  Rates  From  a  Step  Aileron  Input 


Figure  5.12:  Measured  and  True  Angular  Rates  From  a  Step  Elevator 
Input 
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Figure  5.13:  Measured  and  True  Angular  Rates  From  a  Step  Rudder  Input 
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C.  PITCH,  ROLL,  AND  HEADING  SENSOR  MODELING 

The  last  group  of  sensors  to  be  modeled  are  those  which  measure  the  pitch,  roll, 
and  heading  angles.  The  pitch  and  roll  angle  measurements  are  made  with  liquid 
pendulous  devices.  These  are  devices  that  have  an  electrolyte  contained  in  a  vial, 
and  by  measuring  the  capacitance  changes  of  the  vial  as  the  electrolyte  moves  in 
response  to  aircraft  angle  changes,  the  pitch  and  roll  angles  can  be  determined.  The 
heading  sensor  is  a  magnetic  heading  sensor. 

The  primary  concern  with  angular  measurements  is  that  accurate  readings  are 
obtained,  even  when  the  aircraft  is  experiencing  linear  and  angular  accelerations. 
When  these  errors  cannot  be  corrected,  the  control  system  must  be  able  to  compen¬ 
sate  for  the  errors.  An  inclinometer  is  typically  modeled  as  a  pendulum,  shown  in 
Figure  5.14  attached  to  a  block,  which  can  be  considered  to  represent  an  aircraft. 
The  equations  describing  the  motion  of  the  pendulum  are  derived  in  detail  later  in 
this  section.  The  transfer  function  for  the  pendulum  inclinometer  can  be  represented 
as  a  second  order  transfer  function  with  =  0.8  Hz  and  (  =  0.5.  [Ref.  WAT  93] 

0  5.032 

0  ~  s2  + 50.3s +  5.032' 

1.  Error  Modeling 

The  inclinometers  are  subject  to  several  sources  of  errors,  from  both  linear 
and  angular  accelerations,  and  from  mechanical  imperfections.  Mechanical  errors  are 
listed  in  Table  5.3,  and  modeled  as  shown  in  Figure  5.15.  Errors  due  to  angular 
velocity  can  be  compensated  for  with  a  complementary  filter.  Selecting  the  time  con¬ 
stants  appropriately  will  allow  direct  measurements  of  the  angle  from  the  inclinometer 
at  low  frequencies,  and  from  integrated  angular  rates  at  higher  frequencies  where  the 
inclinometer  is  inaccurate. 


70 


Figure  5.14:  Simple  Pendulum  Inclinometer 


TABLE  5.3:  INCLINOMETER  AND  HEADING  SENSOR  CHARAC¬ 
TERISTICS 


Pitch  and  Roll  Range 

±50  c leg 

Pitch  and  Roll  Bandwidth 

1/2  Hz 

Pitch  and  Roll  Accuracy 

0.2  deg 

Heading  Range 

±180  deg 

Heading  Accuracy 

3.0  d(g 

Heading  Repeatability 

0.5  deg 

Heading  Linearity 

0.5% 

It  was  determined  that  the  effects  of  linear  acceleration  were  also  important 
to  model,  in  order  to  model  these  effects  it  was  necessary  to  derive  a  transfer  function 
from  th'  aircraft’s  linear  acceleration  to  the  angie  caused  by  that  acceleration.  The 
starting  point  was  to  model  the  inclinometer  as  shown  in  Figure  5.14.  The  derivation 
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Figure  5.15:  Inclinometer  Error  Modeling 


results  in  two  coupled  equations,  the  first  of  which  will  be  ignored,  since  the  equations 
of  motion  for  the  aircraft  are  known  and  the  motion  of  the  inclinometer  should  have 
little  effect  on  the  aircraft  motion.  The  second  equation  is  for  the  motion  of  the 
pendulum  as  influenced  by  the  aircraft's  acceleration  and  is  what  was  used  for  the 
linear  acceleration  errors. 

First  it  was  necessary  to  define  the  coordinates  used  to  describe  the  incli¬ 
nometer  system  as  q  =  (x.0)  (see  Figure  5.14).  Next,  since  Lagrangian  methods  were 
used  for  the  derivation,  the  kinetic  energy,  T,  of  the  system  was  defined  as 

T  =  1  /2MP  +  1  /2m (x  +  /(?cos0)2  -  !  /'2m(l0  sin#)2 

=  1/2MP  +  1/2  +  1/2  mP2y,  (5.6) 

where  PT  and  Py  represent  the  position  of  the  pendulum  in  Cartesian  coordinates. 


V 


The  potential  energy,  V,  of  the  pendulum  ran  be  written  as 

V  -  V"m  —  mg l (l  —  cos#), 


(5.7) 


where  M,  m,  and  1  are  the  mass  of  the  aircraft,  mass  of  the  pendulum,  and  distance 
from  the  aircraft  c.g.  to  the  pendulous  mass.  The  terms, P,  and  # ,  represent  the 
position  vector  to  the  pendulous  mass  and  the  angle  made  by  the  pendulum  with  the 
vertical  plane.  The  governing  equation  for  Lagrangian  dynamics  is 


£dC  _ 

dt  dq,  dq , 


(5.8) 


where  C  is  the  Lagrangian  operator,  C  =  T  —  V\  and  Q,  represents  the  non- 
conservative  forces  acting  on  the  body.  After  some  rearranging. 


C  =  1/2  (M  +  m)i2  +  mlxQcosO  +  1/2  m/2#2  —  mgl(  1  —  cos  9)  (5.9) 


The  result  in  Equation  5.9  can  be  substituted  into  Equation  5.8  resulting  in 


dC 

dx 

dC 

de 

dC 

dx 

dC 

do 

d_dC 

dt  dx 

d_dC 
dt  d6 


=  (M  +  m)x  +  mlOcosO 
=  mlxcosO  +  PO 
=  0 

=  —mlxO  sin#  —  mgl sin# 

=  (M  +  m)x  +  m/cos##  -  ml  sin##2 
=  ml  cosOx  —  m/ sin#  #  +  inl20. 


When  these  partial  derivatives  are  substituted  into  Equation  5.8,  the  resulting  equa¬ 
tions  of  motion  for  the  pendulum  are 

(A/  +  m)jr  +  ml  cos##  -  ml  sin##2  =  Q\ 
ml  cos  Ox  +  ml20  +  mgl  sin#  =  Q2. 
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where  Q\  and  Q2  are  terms  representing  damping  in  the  system.  The  term  Q 1  repre¬ 
sents  damping  on  the  aircraft  for  Q j  =  —  8Ti  and  the  term  Q2  represents  the  viscous 
damping  of  the  pendulous  mass,  Q 2  =  —j3r6.  The  equations  can  be  linearized,  where 

[M  4-  m)x  +  mid  +  /3xx  =  0  (5.10) 

mix  +  ml26  + /3r9  +  mglQ  =  0.  (5.11 ) 


Now  it  is  apparent  that  Equation  5.10  is  completely  determined  by  the  aircraft  equa¬ 
tions  of  motion  that  have  already  been  modeled.  A  Laplace  transform  of  Equa¬ 
tion  5.11  can  be  performed  to  find  the  desired  transfer  function 


0  _  -1// 


(5.12) 


where  the  similarity  to  the  standard  second  order  transfer  function  is  noted,  making 
the  term  g/l  ~  The  term  — 1  //  in  the  numerator  can  be  solved  for  by  substituting 
g  =  32.174  f/s 2  and  u ;n  =  5.027  rad/stc.  or  0.8 Hz,  with  the  result  /  =  1.27  /.  This 
result  is  then  substituted  into  the  block  diagram  for  the  inclinometer  models  shown 
in  Figure  5.15.  Responses  for  step  inputs  are  shown  in  Figures  5.16,  5.17,  and  5.18. 
In  the  lateral  cases,  there  is  very  little  difference  between  the  measured  and  actual 
angles.  The  longitudinal  case  is  quite  different.  In  Figure  5.16.  there  is  a  considerable 
difference  between  the  actual  pitch  angle.  0.  and  the  pitch  angle  measured  with  the 
inclinometer  model.  This  difference  the  linear  acceleration  of  the  aircraft  as  a  step 
elevator  input  is  applied,  and  must  be  compensated  for  before  a  reliable  control  system 
can  be  developed. 


Angles,  Phi  and  Psi,  radians  Angle  Theta,  radans 


Figure  5.16:  Response  to  a  Step  Elevator  Input 


0  5  10  15  20  25  30  35  40  45  SO 

Time,  seconds 

Figure  5.17:  Response  to  a  Step  Rudder  Input 


i  o 


The  error  in  inclinometer  readings  due  to  linear  acceleration  has  been  de¬ 
termined  experimentally  by  the  manufacturer,  and  can  be  represented  as  a  second 
order  transfer  function  from  g's  to  0, 


0  157.082 

g  ~  's2  +  157.08s  +  157.082’ 


(5.13) 


where  the  gain,  A',  was  determined  from  an  input  of  1.2  V  at  10  m V/g,  and  the 
resulting  output  of  0.4  V  at  60  mY’/deg.  Thus  the  gain,  A',  was  17.99. 


D.  MODELING  OF  AN  ARBITRARY  SENSOR  PLACEMENT 

The  actual  sensor  placement  in  either  the  Bluebird  or  the  AROD  aircraft  is  not 
at  the  c.g.,  as  was  assumed  in  the  previous  sections.  In  order  to  model  the  actual  linear 
and  angular  accelerations  at  the  sensors,  regardless  of  the  position  on  the  aircraft, 
new  equations  of  motion  must  be  derived.  These  equations  can  then  be  applied  to 
the  sensor  inputs  to  obtain  a  proper  output  from  the  sensor.  First  the  equations 
for  linear  acceleration  will  be  derived  using  Newton’s  third  law  for  conservation  of 
momentum,  then  the  equations  for  angular  acceleration  will  be  derived  using  Euler's 
law  for  conservation  of  angular  momentum.  These  equations  must  be  expressed  in 
terms  of  {B}  since  all  the  information  measured  by  the  sensors  will  be  determined  in 
the  bod)  coordinate  system. 

1.  Linear  Accelerations 

In  the  derivation  for  equations  of  motion  of  an  arbitrary  point  on  the  air¬ 
craft,  the  coordinate  systems  representing  the  inertial  and  body  coordinate  systems, 
similar  to  those  shown  in  Figure  3.1  are  used.  Supplementing  the  derivation  of  equa¬ 
tions  of  motion  for  an  aircraft,  the  motion  of  the  sensor  location  on  the  aircraft,  given 
by  Pq  is  necessary.  In  the  inertial  coordinate  system,  {U},  the  position  of  Q  can 
be  written  as  v Pq  =  LbRbPq  -f  1  Pbo-  The  velocity  is  first  determined  by  applying 


Coriolis’  theorem  from  Equation  3.24  to  obtain  the  first  derivative.  Here, 

UVQ  =  vPq  =  vPbo  +  Ujt(LBRBPQl  (5.14) 

where 

Vj/bRbPq)  =  Bj/bRBPq)  +  X  (1bRBPq)-  (5.15) 

The  terms  1 '  jt  and  Bj-t  refer  to  derivatives  taken  with  respect  to  the  inertial  and  body 
coordinate  systems,  respectively.  The  resulting  expression  for  velocity  is 


UVQ  =  uvB0  +  ;x  &rbpq), 


(5.16) 


where  B  jt{lB RbPq  )  =  0  since  Q  is  a  point  fixed  on  the  aircraft.  Accelerations  are 
derived  by  applying  Equation  3.24  twice  to  Equation  5.16 

L%  =  Bjtll'VBO  +  %  x  (gfl%))  +B  X  ("Vio+Jx  (5.17) 

The  B  jt(-)  term  is  differentiated,  resulting  in 

s|(‘'vBO  r,x  (ubrbpq))  =  Bjt(uVBo)  +  V*B  x  1bRbPq  +;x  vvq.  (s.is) 

Substituting  the  results  from  Equation  5.18  into  Equation  5.17,  the  result  is 

UVQ  =  BJt(UVBO )  +  2  %  x  L%  +  V^B  x  lbRbPq  +  £  x  (£  x  1bRbPq).  (5.19) 

Tiie  desired  result  in  {B}  is  obtained  by  simply  premultiplying  Equation  5.19  by  LBR 
d 


B 


Vq  =  — (° vbo )  +  2  bub  x  bVq  -f  bujb  x  bPq  +  bu)B  x  {bub  x  bPq),  (5.20) 
at 


where  the  identity  [Ref.  Sp  89] 


v 


R(UB  X  VVQ)  =  (BRux'b)  X  (?RuVq 


B 


(B  Dt'i 


was  used.  This  result  can  i.ow  be  substituted  into  Equation  3.27  and  the  resulting 
expression  can  be  solved  for  jt(Bi'Bo)- 


V 


2.  Angular  Accelerations 

No  further  work  is  needed  to  derive  the  expressions  for  angular  velocity 
and  acceleration  at  a  arbitrary  point,  since  for  a  rigid  body,  these  quantities  remain 
the  same  anywhere  on  the  body. 

ubo  —  uq  VQ  €  {B} 

^>bo  =  uq  V<2  €  {5}, 

since  m  x  u  =  0. 
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VI.  CONCLUSIONS  AND 
RECOMMENDATIONS 

A.  CONCLUSIONS 

Based  on  the  data  presented  in  this  thesis,  the  following  conclusions  are  made: 

•  High-fidelity  models  of  both  the  AROD  and  Bluebird  aircraft  were  implemented 
using  SlMULINK  software.  Use  of  the  block  diagram  structure  of  the  model 
allows  changes  to  be  easily  made  and  requires  no  programming  ability,  other 
than  the  use  of  MaTLAB. 

•  These  models  accurately  represent  the  sensors,  actuators,  and  engines  associated 
with  the  particular  aircraft. 

•  Expected  errors  in  the  accelerometers  and  rate  gyros  are  very  small.  Pitch 
errors  from  the  inclinometers,  due  to  linear  accelerations,  will  be  much  greater. 

•  The  models  that  were  developed  only  represent  one  flight  condition.  The  AROD 
was  modeled  only  in  a  hover  and  the  Bluebird  was  only  modeled  in  a  cruise 
condition.  Further  work  on  the  control,  guidance,  and  navigation  systems  will 
be  concentrated  in  these  flight  regimes.  The  data  files  for  a  given  flight  condition 
are  easily  replaced  with  tables,  when  more  test  data  are  available. 

B.  RECOMMENDATIONS 

Based  on  the  conclusions  presented  above,  and  the  problems  encountered  during 
this  study,  the  following  recommendations  are  made. 

•  The  Department  of  Aeronautics  and  Astronautics  should  update  the  UNIX  labs 
to  include  SlMULINK  with  MaTLAB  4.0.  This  will  allow  increased  instructional 
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use  in  flight  dynamics,  controls,  and  avionics  courses.  Additionally,  a  practical 
research  tool  would  be  available. 

•  The  work  in  this  thesis  must  be  modified  to  include  the  aerodynamic  charac¬ 
teristics  of  the  Archytas  aircraft.  Achieving  this  step  will  require  wind  tunnel 
testing  for  the  Archytas  in  all  the  expected  phases  of  flight,  especially  the  ex¬ 
tremely  non-linear  transition  phase  from  vertical  to  horizontal  flight. 

•  •  Integration  of  the  quaternion  based  rotation  matrix  should  be  completed  in 
order  to  take  advantage  of  the  superior  qualities  of  quaternions  in  computational 
models. 
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APPENDIX  A:  MATHEMATICAL 
PROPERTIES 


In  this  appendix,  some  of  the  mathematical  properties  used  in  the  text  are 
described. 

A.  CROSS  PRODUCT  PROPERTIES 

In  this  section,  important  properties  of  the  cross  product  and  cross  product 
matrices  are  described  The  cross  product  between  two  vectors. 


A  = 


o.x 

1 

'  bT  ' 

ay 

and  B  = 

by 

.  a>  . 

_  bz  _ 

is  defined  by 


Ax  B  ~ 


Qybz  Q  zby 

a:bT  —  axb: 

a  jby  Uy  6j- 


Properties  of  the  cross  product  are: 


A  x  B  =  ( Ax)B .  where  Ax  = 


and  is  called  a  cross  product  matrix 


az 

—au 


0  -a:  a, 

0  -a2 


B 


R{V  x  U)  =  {gRV)  x  {gRU)  if  g R  is  a  rotation  matrix  with  V  and  U  in  the 


same  coordinate  system.  That  this  matrix  distributes  across  the  cross  product 
is  obvious  since  rotation  matrices  preserve  space  geometry. 

•  gR(V  ■  U)  —  (gRV)  ■  (gRU)  for  the  same  reasons  as  above. 

•  A  x  (B  x  C)  =  (A  ■  C)B  -  (A  •  B)C 

82 


•  A  x  (B  x  C)  =  B  x  (A  x  C)  +  C  x  (B  x  A) 

•  A  x  B  =  -B  x  A 

•  —Ax  =  AT  x 

B.  DERIVATIVES  OF  VECTORS 

For  any  free  vector,  bVq  (i.e.  a  velocity  vector  and  any  rotation  matrix,  BR,  the 
derivative  of  the  velocity  of  Q  computed  in  { J5 }  and  expressed  in  {A}  and  denoted 
as  a(bVq)  is  given  as  [Ref.  Sil  91] 

=  2,{abrbvb) 

=  ABRBVQ  +  ABRjt(BVa) 

=  {ABRBR)ABRB\'Q  +  *BRjt(BVq) 

since  as  shown  in  [Ref.  Cra  86], 

4(  VV'o))  =  aQb  x  a(bVq)  +  abR^-(bVq) 

at  ■  at 

so  then 

ABRBAR  =  An  sx 

The  same  process  can  be  carried  out  for  an  angular  velocity  vector  aQb- 

i{As>B>  =  jt[BRB('nB)) 

=  AnBxAn„  +  ABRjt(B(AnB)) 

=  ABRjt(B(AnB)) 

And  if  the  origins  of  the  coordinate  systems  {A}  and  { B }  are  coincident,  the  derivative 
can  be  expressed  as 

B  •  _  ^  ,B  ,  x 

<^'b)- 
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For  an  applied  vector  bPq  (i.e.  the  position  of  point  Q  in  the  {5}  coordinate  system) 
and  a  rotation  matrix  gR,  the  time  derivative  of  the  position  vector  of  Q.  expressed 
in  {A}  ,  jt(A Pq)  as  a  function  of  it’s  derivative  computed  in  {B}  is  given  by  [Ref. 
Sil  91] 

jt(APc)  =  ~[abR  bPq  P  A Pb°) 

=  ABRBPQ  +  ABRjt(BPQ)  +  AVB 
=  AnBx(ABRBPQ)+ABRjt(BPQ)  +  AVB 
Therefore  the  velocity  of  the  point  Q  can  be  expressed  in  {/f)  as 

aVq  =  A!lB  x  (; aR  bPq)  4  abR  bYq  4  aVb 
or  expressing  the  velocity  of  Q  in  the  {B}  coordinate  system 

B(AVQ)  =  B(AnB)xBPQ  +  BVQ  +  B(AVQ) 

In  the  case  where  the  origins  of  {.4}  and  { B }  are  coincident  then  the  resulting  ex¬ 


pression  is 


B  B  ,  Bn  i  B i /  .  J9 

Vq  —  uJg  X  Pq  A-  VQ  +  Vg. 


C.  EQUATIONS  USED  FOR  LINEARIZATION 

For  any  vector  equation,  H(c)  =  a  x  b  where  c  =  [a  b]  and  a,  b,  and  c  are  all 
vectors,  the  Taylor  series  approximation  can  be  w'ritten  as 

3  H 

H(c)  =  H(co)  +  ~\c=CoSc+H.O.T. 
oc 


3H  3H  c  3H  CL 
dc  ~  da  la=a^G+  db  ' h=b°6b 


Now  dlijda  can  be  written  as 


3H  3{a  x  b)  d(  —  b  x  a)  d(—b)  3a 

37  =  —a -  =  “ sT-  =  -gp  x  0  “ 6  x  77  =  -bx 
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Using  fhe  same  relations,  dH/db  can  be  written 


dH  _  d(-b  x  a) 
db  ~  db 


—  ax 
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APPENDIX  B:  NUMERICAL  RESULTS 


A.  AROD  RESULTS 

The  following  results  are  from  the  numerical  linearization  of  the  kinematic  equa¬ 
tions  of  motion  for  the  AROD.  Using  the  trim  command,  based  on  an  initial  O0  =  tt/2, 
resulted  in  the  state  vector  and  input  vector  of: 

0 
0 
0 
0 

x  =  0 

0 
0 

1.5700 

0 

The  linmod  command,  [a.b,c.d]=linmod('baslnmod'.x.u).  produced  the  following  lin¬ 
earized  system: 


'  0 

0 

0 

0 

0 

0  ' 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1. 

5024 

_  0 

0 

0 

0 

1. 

5112 

0  . 

■  1 

0 

0 

0 

0 

0  ■ 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

b 

— 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

.  0 

0 

0 

0 

0 

1  . 

and  the  c  and  d  matrices  were  empty,  since  no  outputs  were  defined. 

The  following  results  are  for  the  numerical  linearization  with  gravitational  ef¬ 
fects  added  to  the  2-3-1  Euler  angle  kinematic  equation  model.  The  trim  com- 


and  u  = 


0 

0 

0 

0 

0 

0 
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mand.  using  the  same  initial  conditions  as  before,  resulted  in  the  following  vectors: 
[x,u]=trim('gravlmod\xO,[],[].ix,  [].[]) 


Linearizing  the  system,  based  on  the  state  and  input  vector  found  by  trimming  at 


the  desired  conditions  resulted  in:  [a.b.c,d]=linmod('gravlmod\x.u) 


0  0.00  0  0  0.00 

-0.00  0  0.00  -0.00  0 

0.00  -0.00  0  -0.00  -0.00 

0  0  0  0  -0.00 

0  0  0  0.00  0 

0  0  0  -0.00  1.5112 

0  0  0  1.00  -0.00 

0  0  0  0  1.00 

0  0  0  0  0.00 


0.00  0  -0.0255  0.0002  ‘ 

0.00  0.0256  -0.00  32.1740 

0  -0.00  -32.1740  0 

-0.00  0  0  0 

1.5204  0  0  0 

0  0  0  0 

0.00  0  0  -0.00 

-0.00  -0.00  0  0.00 

1.00  -0.00  0  0  . 


■  1.00  0  0  0  0  O' 

0  1.00  0  0  0  0 

0  0  1.00  0  0  0 

0  0  0  1.00  0  0 

0  0  0  0  1.00  0 

0  0  0  0  0  1.00  . 


Trimming  the  full  EOM  model  at  hover  conditions  resulted  in  the  following: 
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The  subsequent  linearization  for  the  complete  model  yielded: 


■  0 

0 

0 

0 

0 

0 

0 

0.0002 

0.0002  ■ 

0.00 

0 

0 

0.00 

0 

0.00 

0.00 

-0.00 

32.1740 

0 

0 

0 

0.00 

0 

0 

-0.00 

-32.1740 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.00 

0 

-1.5174 

0 

0 

0 

0 

0 

0 

-0.00 

1.5082 

0 

0 

0 

0 

0 

0 

0 

1.00 

-0.00 

0 

-0.00 

0 

0 

0 

0 

0 

0 

1.00 

-0.00 

0.00 

0 

0.00 

0 

0 

0 

0 

0 

1.00 

0.00 

0 

0  . 

and 


0 

0 

0 

0.0112 

0 

0 

0 

0 

0 

0 

0 

0 

0 

,*> 

24.8194 

0.0003 

-6.6192 

0 

0 

-0.00 

0 

-6.5791 

0 

-0.00 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Computation  of  the  eigenvalues  gave: 


eig(a) 


0  ' 

0 

0 

0 

0 

-0.00 
0  +  1.5128? 

0  -  1.5128? 

0 


Now  consider  the  quaternion  based  model  Rpm  is  fixed  at  6800  RPM. 


To  hold  the  desired  vertical  attitude,  the  quaternion  states  were  held  fixed: 


it  = 


7 

8 

9 

10 


88 


where  the  state  vector  of  nominal  conditions  was: 


jO  = 


0 

0 

0 

0 

0 

0 

0.7071 

0 

0.7071 

0 


Trim  and  linearization  of  the  quaternion-based  gravitational  model  resulted  in  the 


following  state  and  input  vectors,  as  well  as  linearized  a  and  b  matrices: 


'  32.1740 
-0.00 
-0.00 
0.00 
0.00 
0.00 


a  =  ,  columns  1-6, 


0 

0.00 

0.00 

0 

0.00 

0.00 

0.00 

0 

0.00 

-0.00 

0 

0.00 

0.00 

-0.00 

0 

-0.00 

-0.00 

0 

0 

0 

0 

0 

0.00 

-0.00 

0 

0 

0 

0.00 

0 

-1.6051 

0 

0 

0 

0.00 

1.6062 

0 

0 

0 

0 

0.00 

-0.3536 

-0.00 

0 

0 

0 

0.3536 

-0.00 

0.3536 

0 

0 

0 

0.00 

0.3536 

0.00 

0 

0 

0 

-0.3536 

-0.00 

0.3536 

-0.00 

0.00 

-0.00 

0.00 

-0.00 

0.00 

0.7071 

-0.00 

0.7071 

0.00 
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columns  7-10 


45.5009 

0 

-45.5009 

0 

0.00 

-45.5009 

0.00 

45.5009 

45.5012 

-0.0003 

-45.5012 

0.0003 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.00 

0.00 

-0.00 

0.00 

0 

0.00 

0.00 

-0.00 

-0.00 

0 

0.00 

0.00 

-0.00 

-0.00 

0 

r  i.oo 

0  0 

0 

0 

0  1 

0 

1.00  0 

0 

0 

0 

0 

0  1.00 

0 

0 

0 

0 

0  0 

1.00 

0 

0 

0 

0  0 

0  1.00 

0 

L  0 

0  0 

rv 

U 

0 

1.00  . 

Eigenvalues  are  computed  as: 


eig(a) 


-0.00  +  0.00/ 
-0.00  -  0.00/ 
0.00 
0 

-0.00+  1.6057/ 
— 0.00  —  1 .6057/ 
-0.00 
0.00 

-0.00  +  0.00/ 
-0.00  -  0.00/ 


B.  BLUEBIRD  NUMERICAL  RESULTS 


Results  for  the  numerical  trim  and  linearization  of  the  Bluebird  model  are  given 
below.  The  results  are  from  the  trim  of  the  kinematic  equations  of  motion  for  the 
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Bluebird  model  are  based  on  the  nominal  conditions  of: 

’  72.9954  ' 

0 

6.6757 

0 

tO  =  C 

0 
0 

0.0912 
0  . 

with  the  states,  u,  w,  and  0  held  constant,  it  =  [1  5  8]'.  The  trim  command 

[x,u]=trim(’basic",xO.[],[],ix, [],[])  resulted  in  the  following  state  and  input  vectors: 

'  72.9954  ‘ 

0 

6.6757 
0 

x  =  0  ,  and  u 

0 
0 

0.0912 
0  . 

The  numerical  linearization  of  the  basic  model,  [A,B,C.D]=linmod('basic\x.u),  re¬ 
sulted  in  the  following  A  and  B  matrices: 

'  0  0  0  0  -6.6757 

0  0  0  6.6757  0 

0  0  0  0  72.9954 

A~  0  0  0  0  0 

0  0  0  0  0 

.0  0  0  0  0 

and 

'  1  0  0  0  0  0’ 

0  1  0  0  0  0 

0  0  1  0  0  0 

000100 
0  0  0  0  1  0 

.0  0  0  0  0  1. 

The  C  and  D  matrices  were  empty  since  no  outputs  were  defined. 


0 

-72.9954 

0 

0 

0 

0  . 
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Results  for  trim  and  linearization  of  the  kinematic  mode!  with  gravitational 


effects  added  yielded  the  following  state  and  input  vectors: 


x  = 


’2.9954 

0.00 

6.6763 

0.00 

0.00 

-0.00 

-0.00 

0.0912 

0.00 


.  and  u 


2.8772  ‘ 
0.00 
-31.4605 
-0.00 
0.3945 
-0.00 


The  numerical  linearization  of  the  model  resulted  in  the  following  A  and  B  matrices: 
A  = 


'  -0.0007 

-0.00 

0.0079 

0 

-6.6763 

0.00 

-0.00 

-32.0451 

0 

0.00 

0 

0.00 

6.6763 

0 

-72.9954 

32.0403 

0.00 

0 

-0.0001 

r> 

u 

0.0007 

-0.00 

72.9954 

0 

-0.0002 

-2.8773 

0 

-0.00 

0.00S7 

-0.00 

0 

-0.00 

0.00 

-0.00 

-0.00 

0 

-0.00 

0.00 

0.0005 

0.00 

0 

-0.00 

0.00 

0.0361 

0 

-0.00 

0.0010 

0.00 

-0.00 

-0.00 

0 

-0.00 

-0.00 

0 

0 

0 

0 

1.00 

-0.00 

0.0915 

0.00 

-0.00 

0 

0 

0 

0 

0 

1.00 

0.00 

0.00 

0 

0 

0 

0 

0 

0 

-0.00 

1.0042 

0.00 

-0.00 

0 

and 

'1  0  0  0  0  0' 

0  1  0  0  0  0 

0  0  10  0  0 

a  ~  0  0  0  1  0  0 

0  0  0  0  1  0 

.0  0  0  0  0  1. 

Results  for  trimming  the  entire  model  were  A  = 


'  -0.0622 

0.00 

0.3431 

0 

-1.6187 

0.00 

0.00 

-32.0416 

0 

0.00 

-0.3865 

0.00 

1.7450 

0 

-71.6817 

32.0403 

0.00 

0 

-0.7558 

0.00 

-4.7145 

0.00 

67.1233 

0 

-0.0002 

-2.8771 

0 

0  00 

-0.1455 

0.00 

-5.3700 

0.00 

1.5003 

0.00 

0.00 

0 

0.0155 

0.00 

-0.1907 

0.00 

-3.1266 

0.00 

0.00 

0.0362 

0 

0.00 

0.1418 

0.00 

-1.0589 

0.00 

-0.7970 

0.00 

0.00 

0 

0 

0 

0 

1.00 

0.00 

0.0915 

0.00 

0.00 

0 

0 

0 

0 

0 

1.00 

0.00 

0.00 

0 

0 

0 

0 

0 

0 

0.00 

1.0042 

0.00 

0 

0 
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and 


B  = 


-4.3850 

0 

0 

8.7745 

0.00 

5.6803 

0 

0 

-37.8929 

0 

0 

0 

0.00 

0.6216 

45.9851 

0 

-21.4821 

0 

0.00 

0 

0.00 

-7.1281 

-6.1465 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Complete  model  linearized  with  0O  =  0  instead  of  0O  =  Q0  The  initial  conditions 
are  now: 

73.3000 
0 
0 
0 

j0  =  I  0 
o 
0 
0 
0 

The  state  and  input  vectors  obtained  from  trimming  at  this  state  are 

73.3000 
-0.00 
1.6086 
-0.00 

x  —  |  —0.00  |  and  u  = 

0.00 
0.00 
-0.00 
-0.00 

with  the  linearized  A  and  B  matrices:  A  = 


-0.0181 

-0.00 

0.00 

0.2336 


“  -0.0635 

0.00 

0.3277 

0 

-1.4922 

0 

-0.00 

-32.1740 

0 

-0.00 

-0.3911 

-0.00 

1.6086 

0 

-72.6109 

32.1740 

-0.00 

0 

-0.7572 

-0.00 

-4.7741 

0 

67.9934 

0 

-0.0002 

-0.0002 

0 

0.00 

-0.1471 

-0.00 

-5.4414 

-0.00 

1.5183 

0.00 

-0.00 

0 

0.0151 

-0.00 

-0.1933 

0 

-3.1672 

0 

0.00 

-0.00 

0 

0.00 

0.1440 

0.00 

-1.0578 

0.00 

-0.8114 

0.00 

-0.00 

0 

0 

0 

0 

1.00 

0 

-0.00 

0 

-0.00 

0 

0 

0 

0 

0 

1.00 

-0.00 

-0.00 

0 

0 

0 

0 

0 

0 

0.00 

1.00 

-0.00 

-0.00 

0 

y  3 


-4.5835 

0 

0  , 

8.7745 

0.00 

5.8282 

0 

0 

-38.8681 

0 

0 

0 

-0.00 

0.6252 

47.1717 

0 

-22.0417 

0 

0 

0 

-0.00 

-7.3151 

-6.4345 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

The  eigenvalues  for  the  case  where 


0o  =  0  are 


ag(A) 


0  ' 

-0.5285  +  3.6346/ 
-0.5285  -  3.6346/ 
-5.6291 
—3.9833  +  3.5521/ 
-3.9833  -  3.5521/ 
0.0420 
-0.0191  +  0.4963/ 
-0.0191  -  0.4963/ 


C.  CESSNA  172  RESULTS 


The  Cessnal  172  model  was  trimmed  and  linearized  using  the  state  vector 


x0  = 


i  j 
0 
0 
0 
0 
0 
0 
0 
0 


as  specified  in  [Ref.  Ros  79].  The  states  u,  u\  and  0  were  held  constant  making  the 
term  ix  =  [1  3  8]'. 

Only  the  complete  model  was  linearized,  with  the  results  of  the  trim  routine 
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given  as: 

'  219.00 

0.00 
-0.0394 
-0.00 
x  =  -0.00 

-0.00 
-0.00 
-0.00 
0.00 

The  linearized  model  had  the  following  results  for  the  A  and  B  matrices  A  — 


'  -0.0442 

0.00 

0.0848 

0 

0.0382 

0 

0.00 

-32.1740 

0 

0.00 

-0.1620 

-0.00  - 

0.0394 

0  -217.2141 

32.1740 

-0.00 

0 

-0.2916 

-0.00 

-2.1805 

0  212.5399 

0 

-0.0002 

-0.0002 

0 

0.00 

-0.1313 

0.00  -12.4093 

0.00 

2.5342 

-0.00 

-0.00 

0 

0.0024 

-0.00 

-0.1085 

0  - 

6.0778 

0 

0.00 

-0.00 

0 

-0.00 

0.0462 

0.00  - 

0.3807 

0.00 

- 1 .2600 

0.00 

-0.00 

0 

0 

0 

0 

1.00 

0 

-0.00 

0 

-0.00 

0 

0 

0 

0 

0 

1.00 

-0.00 

0.00 

0 

0 

0 

0 

0 

0 

0.00 

1.00 

-0.00 

-0.00 

0 

and 

'  -6.2509 

0 

0 

3.2255  ' 

-0.00 

19.4571 

0 

0 

-44.3392 

0 

0 

0 

0.00 

4.7446 

57.4954 

0 

B  = 

-39.4919 

0 

0 

0 

-0.00 

-10.2288 

-8.2562 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  . 

and  the  following  eigenvalues 

- 

0  | 

-12.4309 

—0.6947  +  3.3080; 

-0.6947 

-  3.3080/ 

eig{A)  = 

-4.1303  +  4.3895/ 

-4.1303 

-  4.3895/ 

-0.0109 

-0.0209  +  0.1794/ 

L  -0.0209 

-  0.1794/  J 

and  u  = 


0.0001 

0.00 

0.00 

1.0013 
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APPENDIX  C:  PROGRAM  LISTINGS 

A.  AROD  Matlab  ROUTINES 
1.  Main  Routine 

function  accel=mam(vstate,m,rho,A ,R) 

'/,  Function  will  compute  the  accelerations  due  to  the 

'/,  gravitational  forces,  aerodynamic  forces  and  moments, 

*/.  and  control  forces  and  moments. 

'/,  The  values  for  S  , rho  ,m,b  ,and  c  are  used  as  input 

'/,  arguments  to  the  function  call,  and  are  loaded 

X  from  the  workspace.  There  should  be  a  file, 

'/.  loaddata.m  loaded  prior  to  running  the 

'/,  simulation. 

’/.  define  states  in  terms  of  the  input  vector 

u=vstate(l)  ; 

v=vstate(2) ; 

w=vstate(3) ; 

p=vstate(4) ; 

q=vstate(5) ; 

r=vstate(6) ; 

phi*vstate(7) ; 

theta=vstate(8) ; 

psi=vstate(9) ; 
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*/.  Define  the  control  inputs 
de=vstate(10) ; 
dr=vstate(ll) ; 
da=vstate(12) ; 
drpm=vstate(l3) ; 


'/,  this  subroutine  computes  linear  and  angular  accelerations 
7.  given  angular  and  linear  velocities; 

7,  the  input  is  6x1  vector  =  [u  v  w  p  q  r]  ’ 

7.  the  output  is  feedback  part  of  d/dt  [u  v  w  p  q  r]  ' 

7.  the  output  also  includes  the  Euler  angle  derivatives,  based  on 
7.  a  2-3-1  transformation,  for  Ldot ,  used  in  the  function  grav.m. 

v  *  vstate(l :3) ; 
omega  =  vstate(4:6); 
vdot  =  -crpr (omega, v) ; 

[Ib,Ir]  =  inertia; 

7.  compute  the  angular  momentum  due  to  the  body's  inertia,  lb 
Lb  =  lb  *omega; 

7.  compute  the  angular  momentum  due  the  spinning  prop's  inertia,  Ir 
0megaR-drpm*2*pi/6Q ;  7.  angular  velocity  of  the  prop,  rad/sec 

Lr=Ir*[0megaR;0;0] ; 
temp=Lr+Lb; 

omdot  =  -Ib\crpr(omega,temp) ; 
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vstatedot=[vdot jomdot] ; 


'/.  Use  the  Euler  angle  propagation  equations  for  a  2-3-1  rotation  sequence 
Ldot=[p/cos(psi)-cos(phi)*sin(psi)/co..(psi)*q+sin(phi)*sin(psi)/cos(psi)*r 
cos(phi)/cos(psi)*q-sin(phi)/cos(psi)*r; 
sin(phi)*q+cos(phi)*r] ; 

*/.  Given  the  vector  containing  the  state  derivatives, 

*/.  The  function  will  compute  the  forces  and  moments  due  to 

'/.  the  control  derivatives,  Cfd,  where  this 

•/.  is  dCf/dd. 

'/.  The  values  for  rho,A,R,m  are  used  as  input 

'/.  arguments  to  the  function  call,  and  are  loaded 

'/.  from  the  workspace.  There  should  be  a  file, 

'/.  arod.mat  loaded  prior  to  running  the 

'/.  simulation. 

*/. 

'/.  hover  case  V=0,  dimensionalize  the  control  derivatives  based  on 
V,  induced  velocity  through  the  rotor  disk,  Vi  =  sqrt (T/2/A/rho) 

'/.  Calculate  the  quantities  needed  for  the  force 
*/.  calculation. 

*/,  Call  a  function  to  return  the  stability  derivatives  wrt  to  moments 

7.  Could  put  a  call  to  a  lookup  table  here 

*/.  Syntax— Z  =  TABLE2(TAB ,X0 ,Y0)  or  Y  =  TABLE1  (TAB ,X0) 
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Cfd=getcfd; 

*/.  Define  the  derivatives 
Clda=Cfd(4,3) ; 

Cmde=Cfd(5 , 1) ; 

Cndr=Cfd(6,2) ; 

'/,  Calculate  the  Force  due  to  Cfd  derivatives 

X  No  Aerodynamic  Forces  in  a  Hover 

D-0; 

Y=0 ; 

L=0; 


*/.  calculate  the  force  due  to  thrust  in  body  coordinates 
X  USING  THRUST  VS.  RPM  FROM  BOB  STONEY’S  TEST  RUNS 

T=0 . 0297’t‘drpm- 104 . 7 ; 

Vi=sqrt (T/2/A/rho) 

qbar= . 5*rho*Vi"2 ;  X  Vi  is  induced  velocity,  not  forward  speed 
Fout=LD;Y ;L] ; 

Fout=(Fout+[T;0;0] )/m; 

X  Calculate  the  Moment  due  to  Cfd  derivatives 

X  ltr  relates  the  duct  swirl  to  the  moment,  1,  produced  by  thrust 

ltr=-0 . 0542*T-  0.9138; 

l=qbar*A*R* (Clda*da) +ltr ; 

m»qbar*A*R* (Cmue*de) ; 

n=qb^r*A*R*(Cndr*dr) ; 
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Nout=Ib\ [1 ;m;n] ; 

FNcf x= [Fout ;Nout] ; 

*/,  Given  the  vector  containing  the  state  derivatives, 

'/,  and  euler  angles,  the  function  will  compute 

X  the  forces  due  to  gravity  acting  on  the  aircraft. 

y. 

y. 

*/,  Calculate  gravitational  force,  based  on  a  2-3-1  Euler  angle 
y,  rotation  for  position  of  the  aircraft.  Rotation  matrix  is  Ru2body 
y,  Z  for  {U}  is  positive  down  and  X  for  {B}  is  straight  up. 

FNgrav=32 . 174* [-sin (theta) *cos (psi) ; 

sin(theta)*sm(psi)*cos(phi)  +  sin(phi)*cos(theta) ; 
cos(phi)*cos(theta)-sin(theta)*sin(psi)*sin(phi) ; 

0 ;  0 ;  0]  ; 

'/.  Sum  up  the  normalized  forces  and  moments  to  feed  back  into  the  integrator 

vstatedot* [vstatedot+FNcf x+FNgrav ; Ldot] ; 
accel=vstatedot ; 

2.  Supporting  Subroutines 

function  y  *  crpr(omega,x) 
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'/,  this  subroutine  computes  the 
'/.  y  =  omega  X  x 
p  =  omega(l)  ; 
q  =  omega (2) ; 
r  *  omega (3) ; 
t  =  [0  -r  q; 
r  0  -p; 

-q  p  0]  ; 
y  *  t  *  x  i 


crossproduct  of  omega  and  x: 


function  Fgrav=grav(x) 

l  GRAV  will  compute  the  gravitational  force 
*/.  acting  on  the  body,  in  body  coordinates 

g*x(l) ; 
phi=x(2) ; 
theta=x(3) ; 
psi=x(4) ; 

Fgrav=g* [~sin(theta)*cos(psi) ; 

sin(theta)*sin(psi)*cos (phi) +sin (phi )*cos (theta) ; 
cos (phi) *cos (theta) -sin (theta) *s in (psi)*sin (phi)] ; 


function  [ib,ir]  *  inertia 

X  this  subroutine  creates  inertia  matrices  called  ib 
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*/,  and  ir  for  the  body  and  rotor  inertia,  respectively. 

'/,  Arod  hover  case 

ixx  =  1.2312; 

lyy  =  3.9584; 

izz  =  3.9825; 

ixz  =  0; 

irx=. 00898; 

iry=.0045; 

irz= . 0045 ; 

ib  =  [ixx  0  -ixz;0  iyy  0;-ixz  0  izz]; 
ir  =  [irx  0  0 ; 0  iry  0;  0  0  irz] ; 


3.  Data  and  Initialization  Subroutines 


*/.  load  bluebird  data 

A=3. 14; 

*/,  Area  of  rotor  disk,  ft~2 

Vt=712 . 09 ; 

'/,  Rotor  tip  speed,  rad/sec 

m=2 . 6419 ; 

'/.  Mass,  slugs 

R-i; 

X  Radius  of  Rotor  Blade,  ft 

rho* . 002377 ; 

'/.  Air  density,  slug/ft"3 

Uo=0 ; 

'/,  Nominal  Velocity,  ft/sec 

X  Initial  Euler  Angle 

Orientation,  radians 

Phi*0 ; 

Theta*l . 57 ; 

X  Same  as  Steady  State  alpha 

Psi*0; 

Lo* [Phi ; Theta ;Psi]  ; 

7%  Initial  Conditions  to  determine  the  Aircraft 

*/,  Linear  and  Rotational  Velocity  States 

7,  u,v,w  are  computed  from  Uo,  alpha,  and  beta 

'/,  p,q,r  are  assumed  as  zero. 

alpha=Theta; 

beta=0 ; 

Xo=[Uo; alpha; beta]  ; 

7.  Returns  Initial  Conditions  for  the  main 
7*  integrator  in  the  rigid  body  EOM  block 
iO=init_var(Lo ,Xo) ; 

function  Cfd=getcfd 

*/.  Cfd=getcfd  will  return  values  for 

'/,  The  stability  derivatives  for  D,Y,L,l,ia,and  n 

y,  due  to  the  control  inputs. 

'/,  format  for  data  is; 

*/.  [CDde  CDdr  CDda 
y.  CYde  .  .  . 

7  CLde  . . . 

7.  Clde  Cldr  Clda 
7%  Cmde  .  .  . 

X  Cnde  . . .  Cnda 

X  Data  is  non-dimensionalized  by  using  induced  velocity 
X  Vi*sqrt(T/2/A/rho) 
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'/,  for  arod  hover  case 


Cfd=[  0  0  0; 

0  0  0; 

0  0  0 

0  0  1.438; 

-1.233  0  0; 

0  -1.233  0] ; 


function  iO=init_var (Lo ,Xo) 

*/.  INIT_VAR(X)  will  initialize  the  integrators 
7.  initial  states,  iO,  given  the  initial 
7.  conditions  desired. 

7.  Required  initial  conditions  are  the  Euler 
7.  angle  orientation,  total  velocity,  Uo,  initial 
7.  A0A,  and  sideslip  angle,  beta. 

7.  Lo=  [phi  ;  theta ;psi]  ’ 

7.  Xo=  [Uo ;  alpha; beta]  ’ 

All  body  rotation  rates  are  assumed  tc  be  zero 


*/.  Initialize  the  states  u,v,w,p,q,r 
Uo*Xo(l)  ; 
alpha=Xo(2) ; 
beta*Xo(3) ; 
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Ca=cos (alpha) ; 

Sa=sin (alpha) ; 

CB*cos(beta) ; 

SB=sin(beta) ; 

Rwb=[Ca*CB  -Ca*SB  -Sa;SB  CB  0;Sa*CB  -Sa*SB  Ca] ; 
v0=Rwb*[Uo;0;0] ; 
wO= [0 ; 0 ; 0] ; 
i0=[v0;w0;Lo] ; 


function  Qo=initq(lambda) 

'/.  Function  initQ  will  return  values  for 
'/.  the  quaternion  DCM  based  on  a  given 
7.  set  of  Euler  angles. 

7,  Set  for  a  Euler  3-2-1  rotation 

7.  Q(1)=B0 

7.  Q(2)=B1 

7.  Q(3)=B2 

7.  Q(4)=B3 

phi*lambda(l) ; 

theta=lambda(2) ; 

psi*lambda(3) ; 

Qo(l)=.5*sqrt(l+cos(psi)*cos(theta)+sin(psi)*sin(theta)*sin(phi) . . . 
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+cos (ps i ) *cos (phi )+cos (theta) *cos (phi)) ; 

Qo(2)  =  l/4/Qo(l)*(cos(theta)*sin(phi)-sm(psi)*sin(theta)*cos(phi) 
+cos (psi) *sin(phi) ) ; 

Qo(3)=l/4/Qo(l)*(cos(psi)*s  in  (theta)  *cos(phi)+sm(psi)*sin(phi)  .  . 
+sin(theta) ) ; 

Qo(4)  =  l/4/Qo(l)*(sin(psi)*cos(theta)-cos(psi)*sm(theta)*sin(phi) 
+sin(psi) *cos(phi) ) ; 

Qo=Qo  ; 


B.  BLUEBIRD  MatLab  ROUTINES 
1 .  Main  Routine 

function  accel=main(vstate,rho ,b ,c ,S ,m, Xo) 

*/,  Function  will  compute  the  accelerations  due  to  the 

'/.  gravitational  forces,  aerodynamic  forces  and  moments, 

X  and  control  forces  and  moments. 

'/,  The  values  for  S,rho,m,b,and  c  are  used  as  input 

*/,  arguments  to  the  function  call,  and  are  loaded 

X  from  the  workspace.  There  should  be  a  file, 

'/,  loaddata.m  loaded  prior  to  running  the 

*/.  simulation. 

X  define  states  in  terms  of  the  input  vector 
u=vstate(l) ; 
v*vstate(2) ; 
w*vstate(3) ; 
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p*vstate(4) , 
q=vstate(5) ; 
r=vstate(6) ; 
phi=vstate(7) ; 
theta*vstate(8) ; 
psi=vstate(9) ; 

X  Define  the  control  inputs 
de=vstate(lO) ; 
dr=vstate(ll) ; 
da=vstate(12) ; 
dt=vstate(l3) ; 

X  calculate  the  aerodynamic  terms 
Vt=sqrt(u~2+v'2+w~2) ; 
qbar* . 5*rho* (Vt) ~2; 

Ib=inertia; 

X  wind  to  body  transformation 
Rwb=rw2b(vstate,Vt) ; 


X 

X 

X 
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*/,  CHI  will  compute  the  left  hand  side  of  the  state 
'/,  equation.  This  is  the  term  dependant  on  dCf/dxdot. 

'/,  Now  calculate  the  S  matrix  that  non-dimensionalizes 
V.  the  moments.  Also  includes  the  correction  for  Lift  and  Drag 
'/,  acting  in  the  direction  opposite  to  the  positive  coordinate 
direction. 

'/,  Get  the  stability  derivatives  for  forces  and  moments 
Cf xdot=getcf xd ; 

CLad=Cf xdot (3) ; 

Cmad=Cf xdot (5) ; 

%  Twb  is  a  intermediate  step 

Twb= [eye(3) /m  zeros(3) ;zeros(3)  inv(lb)] * [Rwb  zeros(3) ;zeros(3)  Rwb] ; 

'/,  calculate  the  M2  matrix  to  allow  use  of  the  states,  rather  than 
'/.  the  normalized  states.  Only  accounts  for  the  alpha  dot  variables 
*/,  since  the  beta  dot  terms  are  not  ordinarily  computed. 

y,M2*  [0  0  1  0  0  0]*c/2/Vt“2; 
y,Sprime=diag( [-S  S  -S  S*b  S*c  S*b]); 

1.  To  save  some  math  here,  the  product  of  Sprime,  Cfxdot,  and  M2  is: 
PROD® [0  0  0  0  0  0 ; 0  0  0  0  0  0 ; 0  0  -S*CLad*c/2/Vt*2  000; 

0  0  0  0  0  0 ; 0  0  S*c*Cmad*c/2/Vt‘2  0  0  0;0  0  0  0  0  0] ; 

'/,  Calculate  chi 
chil*(eye(6)-Twb*qbar*PR0D) ; 


*/.  Given  the  vector  containing  the  state  derivatives, 

7.  and  euler  angles,  the  function  will  compute 

'/.  the  forces  due  to  gravity  acting  on  the  aircraft. 

V,  Calculate  gravitational  force,  based  on  a  3-2-1  Euler  angle 

*/.  rotation  for  position  of  the  aircraft.  Rotation  matrix  is  Ru2body 

Fgrav=32 . 174* [-sm(theta)  ; 

cos(theta)*sin(phi) ; 
cos(theta)*cos(phi)] ; 


*t  Premultiply  by  the  Chi~-1  term  from  the  left  hand  side 
FNgrav=chil\[Fgrav;0;0;0]  ; 


'/,  Cfx(u)  Given  the  vector  containing  the  state  derivatives, 

*/.  The  function  will  compute  the  forces  and  moments  due  to 

'/.  the  stability  derivatives,  Cfx',  where  this 

%  is  dCf/dx. 

y. 
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7,  Call  a  function  to  return  the  stability  derivatives  wrt  to  moments 

l  Could  put  a  call  to  a  lookup  table  here 

*/.  Syntax— Z  =  TABLE2  (TAB  ,X0  ,Y0)  or  Y  =  TABLE1  (TAB  ,X0) 


Cfx=getcfx; 

CDu=Cf x(l , 1) ; 

CDa=Cf x( 1 , 3) ; 

CYb=Cf x (2 , 2)  ; 

CYp=Cfx(2 ,4) ;  CYr=Cf x (2 , 6) ; 

CLu=Cf x(3 , 1) ; 

CLa=Cf x(3 ,3) ;  CLq=Cfx(3,5) ; 

Clb=Cf x(4 ,2) ; 

Clp=Cf x(4 ,4) ;  Clr=Cfx(4,6) ; 

Cmu=Cfx(5 , 1)  ; 

Cma=Cf x(5 ,3) ;  Cmq=Cfx(5 ,5) ; 

Cnb=Cf x(6 ,2)  ; 

Cnp=Cf x(6 ,4) ;  Cnr=Cf x(6 ,6) ; 

ss=getcf 0 ; 

CD0=ss(l)  ; 

CL0=ss (3)  ; 

Cm0*ss(5)  ; 

Cfd=getcfd; 

'/.  Define  the 

derivatives 

CDde=Cf d( 1 , 1) 

CYdr*Cfd(2,2) 

;  CYda=Cfd(2,3) ; 

CLde=Cfd(3, 1) 

Cldr*Cfd(4,2) 

;  Clda=Cf d(4 ,3) ; 

Cmde*Cfd(5, 1) 

Cndr*Cfd(6,2) 

;  Cnda=Cf d(6,3) ; 
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*/.  Calculate  the  Force  due  to  Cfx'  derivatives 
'/.  And  the  control  derivatives 
7.  in  the  wind  coordinates 

D*-S*qbar/m* (CDO+CDa*w/Vt+CDde*de) ; 

Y*S*qbar/m* (CYb*v/Vt+CYr*r*b/2/Vt+CYdr*dr+CYda*da) ; 

L=-S*qbar/m* (CL0+CLa*w/Vt+CLq*q*c/2/Vt+CLde*de) ; 

7,  calculate  the  force  due  to  thrust  in  body  coordinates 
7.  THRUST  IS  ESTIMATED,  BASED  ON  4.0  HP,  PROP  EFF1CENCY= .  65 

T=15 ; 

Xt=T/m*dt ; 

Fout=Rwb* ( [D ; Y ;L] ) ; 

Fout=Fout+ [Xt ; 0 ; 0]  ; 

7.  Calculate  the  Moment  due  to  Cfx'  derivatives 
7.  And  the  control  derivatives 

l*qbar*S*b*(Clb/Vt*v+Clp*b/2/Vt*p+Clr*b/2/Vt*r+Cldr*dr+Clda*da) ; 

m*qbar*S*c* (Cra0+Cma/Vt*w+Cmq/Vt*c/2*q+Cmde*de) ; 

n*qbar*S*b* (Cnb/Vt*v+Cnp*b/2/Vt*p+Cnr*b/2/Vt*r+Cndr*dr+Cnda*da) ; 

Nout=Ib\(Rwb*[l;m;n] ) ; 
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X  Premultiply  by  the  Chi"-1  term  from  the  left  hand  side 
FNcfx=chil\ [Fout ;Nout] ; 

X 

X - - - 


X  this  subroutine  computes  linear  and  angular  accelerations 
X  given  angular  and  linear  velocities; 

X  the  input  is  6x1  vector  =  [u  v  w  p  q  r]  ' 

X  the  output  is  feedback  part  of  d/dt  [u  v  w  p  q  r] ’ 

X  the  output  also  includes  the  Euler  angle  derivatives 
X  Ldot,  used  in  the  function  grav.m. 

v  *  vstated  :3) ; 

omega  =  vstate(4:6); 

vdot  *  -crpr (omega, v) ; 

temp  =  Ib*omega; 

omdoi  ■  -Ib\crpr(omega,temp) ; 

vstatedot=chil\[vdot ;omdot] ; 

X  Use  the  Euler  angle  propagation  equations 
Ldot= [p+sin(phi) *tan(thet  a) *q+cos (phi) *t an (theta) *r ; 
cos(phi)*q-sin(phi)*r ; 

sin(phi)/cos(theta)*q+cos(phi)/cos(theta)*r] ; 
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vstatedot= [vstatedot+FNcf x+FNgrav ; Ldot] ; 
accel=vstatedot ; 

2.  Supporting  Subroutines 
function  y  =  crpr (omega, x) 

this  subroutine  computes  the  crossproduct  of  omega  and  x 
'/,  y  =  omega  X  x 
p  *  omega(i) ; 
q  =  omega (2) ; 
r  =  omega(3) ; 
t  *  [0  -r  q; 
r  0  -p; 

-q  p  0] ; 
y  -  t  *  x; 

function  Fgrav=grav(x) 

Z  GRAV  will  compute  the  gravitational  force 
7,  acting  on  the  body,  in  body  coordinates 

g=*x(l) ; 
phi*x(2) ; 
theta=x(3) ; 
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Fgrav=g* [-sin(theta) ; 


cos(theta)*sin(phi) ; 
cos(theta)*cos(phi)] ; 


function  lb  =  inertia 

*/.  this  subroutine  creates  inertia  matrix  called  lb 
*/,  for  the  Bluebird  test  aircraft. 

Y.  All  units  are  slug-ft~2 
Ix=10 . 0 ; 

Iy=16 . 12 ; 

Iz=7.97; 

Ib«[Ix  0  0 ; 0  Iy  0;0  0  Iz]  ; 


function  [Rwb] =Rw2b(x , Vt) 

Y,  RWIND2B0DY  Rotation  matrix  for  wind  to  body  coordinate  transformations. 

Ca=x(l)/Vt; 

Sa=x(3)/Vt; 

CB*x(l)/Vt ; 

SB=x(2)/Vt; 

Rwb= [Ca*CB  -Ca*SB  -Sa;SB  CB  0;Sa*CB  -Sa*SB  Ca] ; 
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3.  Data  and  Initialization  Subroutines 


7,  load  bluebird  data 


S=22.38;  */. 

Uo=73.3;  */. 

m= 1.7095;  */. 

b=12 .42;  */. 

c=  1.802;  */. 

rho® .002377 ;  */. 


Planform  Area,  ft'2 
Nominal  Velocity,  ft/sec 
Mass,  slugs 
Span,  ft 

Average  Chord,  ft 
Air  density,  slug/ft~3 


V,  Initial  Euler  Angle  Orientation,  radians 
Phi=0 ; 

Theta®. 0912;  '/.  Same  as  Steady  State  alpha 

Psi®0; 

Lo=  [Phi ;Theta;Psi] ; 


*/.  Initial  Conditions  to  determine  the  Aircraft 

*/,  Linear  and  Rotational  Velocity  States 

V,  u,v,w  are  computed  from  Uo,  alpha,  and  beta 

*/,  p,q,r  are  assumed  as  zero. 

alpha=Theta; 

beta=0 ; 

Xo® [Uo ; alpha ; beta] ; 

y.  Returns  Initial  Conditions  for  the  main 
y,  integrator  in  the  rigid  body  E0M  block 
i0®init_var(Lo ,Xo) 
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function  iO=init_var(Lo,Xo) 

X  INIT_VAR(X)  will  initialize  the  integrators 
X  initial  states,  iO,  given  the  initial 
'/,  conditions  desired. 

'/,  Required  initial  conditions  are  the  Euler 
'/.  angle  orientation,  total  velocity,  Uo ,  initial 
X  AOA,  and  sideslip  angle,  beta. 

X  Lo=[phi ; theta ;psi] ’ 

X  Xo= [Uo ; alpha; beta] ’ 

X  All  body  rotation  rates  are  assumed  to  be  zero 

X  Initialize  the  Euler  angle  DCM 

X  Initialize  the  states  u,v,w,p,q,r 
Uo=Xo(l) ; 
alpha=Xo(2) ; 
beta=Xo(3) ; 

Ca=cos (alpha) ; 

Sa=sin(alpha) ; 

CB*cos(beta) ; 

SB*sin(beta) ; 

Rwb*[Ca*CB  -Ca*SB  -Sa;SB  CB  0;Sa*CB  -Sa*SB  Ca] ; 
v0*Rwb*[Uo;0;0] ; 
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w0=[0;0;0] ; 
iO=[vO;wO;Lo] ; 

function  CfO=getcfO 

X  CfO*getcfO  will  return  values  for 

'/,  the  nominal  values  for  coefficients 

'/.  format  of  input  is  [CDO  CYO  CLO  CIO  CmO  CnO]  ’  ; 

Cf 0= [0 . 03  00.3000]'; 

function  Cfd=getcfd 

*/.  Cfd=getCfd_F(n)  will  return  values  for 
*/.  The  stability  derivatives  for  D,Y,and  L 
'/,  due  to  the  control  inputs. 

7,  format  for  data  is; 

X  [CDde  CDdr  CDda 
X  CYde  . . . 

X  CLde  . . . 

X  Clde  Cldr  Clda 
y.  Cmde  .  .  . 

X  Cnde  . . .  Cnda] 

X  For  the  test  aircraft  Bluebird 
7,  Derivatives  from  DATC0M 
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Cfd= [  .065  0  0; 

0  .0697  0; 

.472  0  0 

0  .0028  .265; 

-1.41  0  0; 

0  -.0329  -.0347]; 

function  Cfx=getcfx 

'/.  Cfx=getcfx(n)  will  return  values  for 
'/,  The  stability  derivatives  for  D,Y,and  L 
*/,  due  to  the  state  vector. 

'/,  format  of  data: 

’/.[CDu  CDb  CDa  CDp  CDq  CDr; 

•/.  CYu  ... 

*/.  CLu  ... 

*/.  Clu  ... 

’/.  Cmu  .  .  . 

*/.  Cnu  ...  ] 

7,  For  the  test  aircraft  Bluebird 

y.  Derivatives  from  DATC0M 

Cfx® [0  0  .188  00  0; 

0  -0.31  0  00  .0973; 

0  0  4.22  0  3.94  0; 
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0  -.0597  0  -.363  0  .1; 

0  0  -1.163  0  -11.77  0; 

0  .0487  0  -.0481  0  -.0452]; 

function  Cfxdot=getcfxd 

%  Cfxdot=getcfxd(n)  will  return  values  for 

l  The  stability  derivatives  for  D,Y,and  L 

'/,  due  to  the  state  vector.  Beta  dot  terms  are  ignored 

'/.  since  they  are  not  normally  determined. 

*/,  format  is: 

'/,[  CDadot 
l  CYad 
*/.  CLad 
%  Clad 
%  Cmad 
'/.  Cnad  ] 

Cf xdot* [0 ; 0 ; 1 . 32 ; 0 ; -4 . 7 ; 0] ; 
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